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SUMMARY 


The present study documents in detail the effect of initial mean velocity 
field on the stability characteristics of longitudinal vortices. The temporal 
stability of isolated multiple cell vortices is considered. The types of 
vortices studied in this work include single cell as well as two and three 
cell vortices. It is shown that cell multiplicity in the vortex core has 
drastic effects on the stability characteristics. On the basis of numerical 
calculations, it is concluded that the growth rates of instabilities in 
multiple cell vortices are substantially larger (two to three fold increases 
are observed) than those of a single cell vortex. It is also determined that 
there is a substantial increase in the effective range of axial and azimuthal 
wavenumbers where instabilities are present, but of most importance is the 
appearance of a variety of viscous modes of instability. In the case of 
vortices, these latter instabilities which highlight the importance of viscous 
forces have never been reported before. For the case of a two cell vortex, 
these effects are discussed in detail. 

A Chebyshev spectral collocation method with staggered grid was developed 
for application to stability calculations for general three dimensional 
swirling flows. The stability problem is formulated in the primitive variable 
form; hence the use of staggered grid eliminates the need for artificial 
pressure boundary conditions. It is shown that this method is robust and 
produces results very accurately and efficiently. 


While validating the spectral collocation algorithm, two new viscous 
modes of instability for Batchelor's vortex were found. These two forms of 
instabilities include an axisymmetric, n = 0, as well as an asymmetric, n = 1 
mode. Here, n is the azimuthal wavenumber of the disturbance. It is shown 
that both of these modes are consistent with the pathology of contrail photos 
of large trailing line vortices as well as with experimental measurements 
conducted in laboratory wind tunnels. Both of these disturbances are 
discussed in some detail. 
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LIST OF SYMBOLS 


a constant for coordinate mapping (Eq. (4.4)) 
constant non-dimensional pressure gradient 
center mode instabilities 
perturbation growth rate (u^/a) 
perturbation phase speed (u> r /a) 
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the radial position where far field boundary conditions 
(at infinity) are applied 

a local radius 

dimensional radial coordinate 

seal ing factors obtained from coordinate mapping 

mean radial velocity 

radial perturbation velocity 

combined mean and perturbation radial velocity 

mean azimuthal velocity 

azimuthal perturbation velocity 

total velocity vector 

perturbation velocity vector 

combined mean and perturbation azimuthal velocity 

mean axial velocity 

axial perturbation velocity 

wall mode instabilities 

a function representing the radially dependent part of the 
mean axial velocity 

combined mean and perturbation axial velocity 
centerline axial velocity 

two constants in the axial velocity of a Batchelor vortex 
(Eq. (4.2)) 

a vector representing the eigenfunctions 
axial coordinate 


Greek Symbol s 


axial wavenumber 
compressibility factor 
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angle as defined in Eq. (3.3) 
flow circulation defined as rV 

a parameter identifying the position of the critical layer 
as defined in Eq. (1.4) 

a similarity parameter for the mean velocity profile 

eigenvalue of the discretized system 

kinematic viscosity 

mean pressure 

azimuthal coordinate 

density 

angular velocity of the flow 
complex angular perturbation frequency 
imaginary part of w 
real part of w 

independent variable in Chebyshev space 

a parameter representing non-parallel effects in mean 
velocity profile 
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Chapter 1 


INTRODUCTION 

The study of the stability of streamwise vortices is an important and 
active area of major current research. Control and manipulation of organized 
longitudinal vortical motions is of interest for several applications such as 
the wake-vortex hazard problem, three-dimensional flow separation control, and 
combustor optimization. The breakdown of leading edge vortices on delta 
wings, which severely reduces lift, and the very stable wing tip vortices shed 
from large commercial aircraft, which determine the flight frequency at 
airports, are two classical examples. 

Although the best (and natural) way to attack these problems is through 
theoretical investigations or full numerical simulations, major difficulties 
exist. The theoretical difficulty with stability analyses was identified in 
the classic work of Howard and Gupta [1]*. They determined that for steady, 
inviscid, incompressible three-dimensional vortices subjected to three- 
dimensional disturbances, no general necessary and sufficient condition for 
stability (or instability) could be found. Hence, to obtain stability 

characteristics, each vortex flow must be studied individually. After twenty- 
five years of active research, the finding [1] still holds today which is a 
source of major disappointment. Aside from the prohibitive cost of each 
numerical simulation, computational results lack the needed spatial 


* 


The numbers 


in brackets 


indicate references. 
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resolution. They are also restricted to low Reynolds number flows, and the 
proper imposed outflow boundary conditions are unclear. Therefore, for 
practical purposes, at the present time, one is forced to rely on the same 
tools which have been developed as far back as seventy years ago and were 
applied successfully to studies of the stability of boundary layer flows. One 
such tool is normal mode analysis. Indeed in the present author s opinion, 
our current understanding of vortex stability resembles the state of affairs 
for boundary layer stability in the 1930's. Our total lack of understanding 
of the importance of such effects as 

(a) initial mean velocity profile 

(b) non-parallel mean velocity profiles 

(c) turbulence 

(d) density stratification 

(e) compressibility 

on the stability characteristics of vortices is widely known and a source of 
concern. Meaningful experimental data are rare because vortices are very 
sensitive to probe intrusion and flow visualization is biased by the strong 
body forces acting on any foreign particles. Also, it is very difficult to 
isolate the effects of each of the parameters mentioned above on vortices, 
much less measure them. Even the stability calculations conducted numerically 
have had their share of problems. The complexity stems from the three- 
dimensional nature of the phenomena and hence the necessity of finding three- 
dimensional mean velocity profiles which are solutions of the Nav ier-Stokes 
equations. Most of the stability calculations to date have used Batchelor's 
vortex [2], and in some cases Long's vortex [3, 4] or solid body rotation 
superimposed on Poiseuille flow in a pipe as the basic unperturbed flow. 
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The present study is directed toward understanding the influence of the 
initial mean velocity profile on the stability of an axial vortex. From an 
engineering standpoint, the research goal is to develop vortex control strate- 
gies via alteration of vortex stability characteristics. Vortex stability can 
be altered through modification of the initial three-dimensional mean-velocity 
field of the vortex. Bushnell (private communications) has put forth the idea 
that once the desired stability characteristics are produced by a particular 
velocity profile, the control problem can be inverted in the sense that the 
wing tip geometry can be tailored to produce the requisite velocity profile. 

First, an overview of some of the important theoretical and computational 
studies regarding hydrodynamic stability of vortices is given. This review is 
by no means exhaustive, but the cited papers are those which have been judged 
relevant to the present research. It must be emphasized here that no attempt 
has been made to relate the current work to the question of vortex breakdown 
as described in the open literature. As has been mentioned by Leibovich [5], 
vortex stability and vortex breakdown are two distinct phenomena (problems) 
which may, in some circumstances, be interrelated. Many reviews have been 
written concerning the breakdown of vortices. The reader is referred to the 
works of Hall [6], Hall [7], Leibovich [5, 8] and Escudier [9] for further 
information regarding vortex breakdown. In Chap. 2, the stability problem is 
formulated and the governing equations, along with boundary conditions, are 
derived. The normal mode analysis is discussed briefly and also the type of 
mean velocity profiles employed in the present work are explained fully. 
Chapter 3 contains the numerical algorithm employed in these calculations. It 
is an extension of the spectral collocation technique developed recently by 
Khorrami, Malik and Ash [10]; hereafter referred to as KMA. The method is 
tested for convergence and accuracy against the results of some well known 
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calculations and the results are reported in Chap. 4. Specifically, two cases 
have been considered; Poiseuille flow in a pipe, and Batchelor s vortex. In 
the latter case, two new viscous modes of instability have been found which 
have not been reported previously. Chapter 5 contains the stability results 
for a variety of mean velocity profiles. The profiles include a single cell 
vortex as well as multiple cell vortices. Finally, the concluding remarks are 
presented in Chap. 6. Here the major results obtained in Chaps. 4 and 5 are 
highlighted and suggestions for future work are given. 

The review in this chapter is broken up into four separate sections. 
Section one contains those stability criteria which are general in nature and 
do not depend on the form of the assumed mean velocity profile. The stability 
results for Hagen-Poiseuil 1 e flow with superimposed solid body rotation is 
discussed in section two. Section three is devoted primarily to the stability 
of Long’s vortex but also includes some other miscellaneous profiles. 
Finally, the stability characteristics of Batchelor's vortex are presented in 
the fourth section. It must be stressed here that, unless otherwise specified 
explicitly, the stability results presented in these sections are for 
inviscid, incompressible fluid. Hence, the viscous results are identified 
explicitly in this review while the inviscid results are not. 

1.1 General Stability Criteria 

The hydrodynamic stability of a vortex was first studied by Lord Rayleigh 
[11]. In his classical paper, he showed (intuitively and using an angular 
momentum argument) that a necessary and sufficient condition for stability of 
a vortex without an axial velocity component was that the square of 
circulation must increase outwards. That is 

iU ^->0 ( 1 . 1 ) 

dr 
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everywhere. Here, r is the circulation and r is the radial coordinate. 
Rayleigh's result was proved in a rigorous manner mathematically by Synge 

[12], who showed that stability of a heterogeneous vortex is assured if and 
only if 


pp (pr 2 ) > 0, 


( 1 . 2 ) 


where p is the density. It is clear from Eq. (1.2) that in the case of a 
constant density fluid, Rayleigh's criterion is obtained. Much later, Howard 
and Gupta [1] were able to give a sufficient condition for the stability of a 
vortex with an axial velocity component (subject to axisymmetric 
disturbances). They found that for stability the local inequality 


1 dr 2 r dW, 2 .. 

“T dT^HF > 1/4 


(1.3) 


must be satisfied everywhere. Here, W is the axial velocity. It is obvious 
from the above criterion that axial shear has a destabilizing effect. It must 
be emphasized that violation of (1.3) does not necessarily mean instability. 


For non-ax i symmetr ic disturbances, Howard and Gupta found a sufficient 
condition for stability in terms of axial and azimuthal wavenumbers. They 
stated that this criterion is always violated for sufficiently small axial 
wavenumbers, and while this does not imply instability, it has led them to 
suggest that no general necessary and sufficient condition is obtainable. As 
indicated earlier, their statement still holds and thus far no new general 
criterion has been obtained. 


1.2 Stability of Hagen-Poiseuil le Flow with 
Superimposed Rigid Rotation 

The stability of rotating pipe flows has been studied by many 
investigators. Using the radial disturbance eigenfunction equation derived by 
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Howard and Gupta, Pedley [13] showed (interestingly in an analytical and 

W 0 

tractable manner) that for a very small Rossby number (e = where q q is 

o o 

the angular velocity of the pipe wall, R 0 is the pipe radius and W Q is the 

centerline axial velocity) the flow is unstable to non-axisymmetric 

disturbances of sufficiently large axial wavelength. He found that, although 

both solid body rotation and Poiseuille flow in a pipe are stable with respect 

to infinitesimal disturbances, the superimposed combined flow (which resembles 

a helical vortex) is highly unstable. These disturbances are helical 

traveling waves which wrap around the mean flow in the opposite direction to 

the streamlines. The growth rate of the most rapidly growing disturbance is 

2e q . In a follow-up paper, Pedley [14] showed that for viscous 

0 

disturbances, the critical Reynolds number, Re c , has a value of 82.9 
corresponding to azimuthal wavenumber, n = -1. The critical Reynolds number 
increases as the value of azimuthal wavenumber , n, increases. The 
disturbances are stationary relative to the rotating frame of reference, and 
as the Reynolds number increases, the wavenumber of the most rapidly growing 
disturbance also increases. Pedley's findings were confirmed by Metcalfe and 
Orszag [15] who employed a spectral tau method; and later on they were 
reconfirmed by KMA. 


Masl owe [16] extended the work of Pedley [13] without making any 

assumption as to the magnitude of the Rossby number, e. He showed that the 
most unstable modes have negative azimuthal wavenumbers, n. They spiral in 
the same direction as the basic flow rotation (note this is in contradiction 
to Pedley's finding which is caused by different interpretations of what 

direction a wave with negative azimuthal wavenumber spirals) but propagate 

upstream in the axial direction with an axial phase speed 0(e _1 ). The 

amplification factor and the axial wavenumber, a, of the fastest growing 


ll 



disturbance peaked at a finite Rossby number e of 0(1) while the growth rate 
showed little variation with n for |n| > 2 at finite values of e . These 
instabilities persist at very large values of e and have larger growth rates 
than suggested previously. Mas! owe was also able to prove that a necessary 
condition for instability could be predicted using a parameter, Y defined by 

Y r = aeQ 0 W + n Q q - ay , (1.4) 

where ay is the real part of the growth rate, u. If Yp was zero anywhere 
in the interval (0, R 0 ), the flow was unstable. Condition (1.4) is a formal 
statement for the existence of a critical layer in inviscid swirling flows. 

The neutral stability curves for rotating Poiseuille flow subjected to 
viscous disturbances is given by MacKrodt [17]. He found that the first mode 
which becomes unstable is for the case of n = - 1 with a critical Reynolds 
number of Re^ - 82.88 which is exactly the value obtained by Pedley. However, 
MacKrodt's major observation was the fact that at high Reynolds numbers, the 
amount of rotation required for destabilization is much smaller than had been 
expected. 

In a very extensive study by Cotton and Sal wen [18], MacKrodt's results 

were extended to much higher Reynolds numbers. Re, and rotation rates, Q . In 

o 

their viscous calculations, they found that the instability region is pushed 
to higher Re and Q q when a > 2. The disturbances are center mode at least 
for high a; that is most of their energy is confined near the centerline of 
the pipe. The bulk of their study was directed toward the n = -1 mode where 
they calculated the neutral stability curves for the primary mode and for 
several higher modes (the existence of these higher modes were first 
established by Duck and Foster [19] for the case of a trailing line vortex; we 
will come back to this point repeatedly as the review proceeds). Cotton and 
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Sal wen found out that although the upper branches of the neutral curves are in 

excellent agreement with the results of Pedley the lower branches behave in a 

more complicated manner. That is, there were many kinks, loops, points of 

degeneracy, branching, and mode crossing or jumping associated with the lower 

branch which had not been identified previously. 

A study of the asymptotic nature of the neutral modes in the case of 

| n | >> i was attempted by Maslowe and Stewartson [20] using perturbation 

theory. Expanding the eigenfunctions in descending powers of n, they showed 

that the neutral disturbances were of the wall mode type and there are 

infinite numbers of neutral stability modes for any n and e, which is 

consistent with the results of Pedley [13] and Duck and Foster [19]. However, 

the neutral wall modes obtained by [20] are not the limiting inviscid modes 

(as Re •*> “>) associated with the viscous modes obtained by Cotton and Salwen. 

This important discrepancy with the viscous results of [18] led [20] to 

suggest that the role of viscosity in swirling flows is more subtle than 

previously thought, and that it alters the shape of the perturbation's 

eigenfunction significantly. For the case of n = -1, Maslowe and Stewartson 

[20] resorted to a numerical method and found the eigenfunction structure 

differed significantly from the modes with large n. At finite values of e, 

the critical layer associated with n = -1 was close to the centerline. 

However, as e + °°, the critical layer moved very slowly toward the pipe 

wall. Finally, Maslowe and Stewartson postulated from their results that if a 

swirling pipe flow meets the conditions 

V(R q ) * 0, (1.5a) 

and 

(R ) * 0, U-5b) 

dr v o' 

then the flow must be unstable. 
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The results for neutrally stable disturbances were obtained in a closed 
form by Ito, Suematsu, and Hayase [21]. However, rather than assuming 
Poiseuille flow, they employed a uniform axial velocity. They found that 
there are an infinite number of higher modes for each distinct set of flow 
parameters. Looking for stationary waves, they found the range of the Rossby 
number for which there exists a stationary axisymmetric mode (n = 0) to be 

0 < E < 0.522 . (1.6) 

However, in the case of non-axisymmetric modes, they stated that no such 
restriction on the value of Rossby number exists. Since they never looked for 
unstable disturbances, no conclusions can be made with regard to the 
possibility that a pipe flow might become unstable even though it does not 
meet the requirement specified in (1.5b). 

1.3 Stability of Long's Vortex 

The similarity solution for a bathtub vortex due to Long [3, 4], has been 
employed by several researchers to obtain the stability characteristics of 
tornado type vortices. Such a vortex was first used by Burggraf and Foster 
[22]. However, they had great difficulties with their numerical approach (a 
local shooting method) in obtaining converged solutions. They attributed the 
failure of the numerical method to the algebraically (as contrasted with 
exponentially) decaying axial velocity W. An algebraically decaying mean flow 
requires a large integration radius at which the far field boundary conditions 
are imposed, which in turn causes great difficulties with the integration step 
size. The dependency of the eigenvalues on this radius will be explained in 
greater detail in Chap. 4. Switching to an exponential velocity profile, 
Burggraf and Foster obtained converged solutions and appear to be the first to 
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report, numerically, the existence of a spectrum of unstable eigenvalues (a 
result which generally is credited to Duck and Foster [19]). They found that 
as in the case of a rotating pipe, the most dangerous disturbances have an 
azimuthal wavenumber, n = -1. 

The calculations of [22] were extended by Foster and Duck [23]. 
Employing a global finite difference technique similar to that of [19], they 
encountered no difficulty in obtaining converged results for Long's vortex. 
They showed that only the negative azimuthal wavenumbers are unstable and the 
most unstable mode is approached in the limit as n + - ®. This result, which 
apparently holds for all inviscid swirling flows, was first reported by 
Lessen, Singh, and Pail let [24] for the case of a trailing line vortex. The 
general belief in the stability community is that viscosity would effect these 
disturbances so that the most unstable mode would occur at a finite value of 
n. Foster and Duck mentioned that every effort to obtain instabilities for n 
> 0 failed; and that a vortex with a single cell core is more unstable than a 
two cell vortex. It is also apparent from their growth rate plots that in the 
case of the Long vortex, inviscid instabilities are an order of magnitude 
weaker than for the case of a trailing line vortex. 

Staley and Gall [25] employed a vortex which has continuous velocity 
components, but whose vorticity field has many discontinuities. Their work is 
mentioned to show that even with such a crude approximation for the mean field 
of a vortex, they obtained results that in general substantiate the findings 
of [19, 23, 24]. That is, the solutions indicate an infinite number of higher 
modes and the most unstable mode is approached as n -*• - ». They also 
obtained an instability criterion (for large n) in terms of the vorticity 
field which appears to be another version of the sufficiency condition first 
reported by Leibovich and Stewartson [26]. 
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1.4 Stability of Batchelor's (Trailing Line) Vortex 

The stability of the mean velocity profile of a trailing line vortex 
(Batchelor [2]) was studied by Lessen et al . [24] with respect to 

infinitesimal, non-axisymmetric disturbances. It was found that negative 
azimuthal wavenumbers are destabilized by the addition of swirl. Lessen et 
al . discovered that the stability characteristics of the vortex are strongly 
dependent on the value of q, a swirl parameter which is related to the ratio 
of the maximum swirl velocity to the maximum axial velocity (this parameter is 
defined in Chap. 4). All wavelengths appear to become damped, and the flow is 
completely stabilized at a value of q slightly greater than 1.5. They found 
that the most unstable wave is obtained in the limit n ->■ - <*>, while the 
axial wavenumber associated with maximum growth rate also increases without 
bound. Following their earlier work. Lessen and Pail let [27] performed a 
viscous stability calculation for the trailing line vortex. They obtained 

results similar to those of Pedley [14] which indicated that the critical 
Reynolds number increases as |n| increases, with the first unstable mode 

corresponding to n = -1. Their calculations of neutral stability curves 

showed that the values of the critical wavelength and critical Reynolds number 
are not very sensitive to the exact value of q, and at large q, stability 
is insured for any wavelength and Reynolds number. However, their results 

are limited to a small range of Reynolds numbers (Re < 150) which later will 
be shown to be inadequate. Duck and Foster [19] extended the work of Lessen 
et al . [24] and found a continuous spectrum of unstable modes. The spectrum 
contains infinite numbers of higher modes for each set of particular values of 
a, n, and q. They showed that the number of unstable modes depends directly 
on the number of grid points, N, so that as N increases the number of these 
unstable modes increases as well. It is shown by KMA that although they are 
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i n v i sc id in nature, these higher modes exist at much lower Reynolds numbers 
than expected. 


The asymptotic behavior of the most unstable mode for large azimuthal 
wavenumber (|n| >> 1) was obtained by Leibovich and Stewartson [26]. They 
showed that as n + - the maximum growth rate, , approaches a constant 
value, u>. . These disturbances are ring modes; that is they are confined to 
the vicinity of a local radius -- say r = r Q . There are also multiple 
critical layers where y p = 0; and it was pointed out by [26] that both the 
ring modes and the wall modes obtained by Maslowe and Stewartson [20] may 
coexist. As will be shown later, this indeed is the case for multiple cell 
vortices where wall modes, ring modes and center modes all coexist. Leibovich 
and Stewartson were able to confirm the results of [19] and showed that as one 
gets nearer to the neutral stability curves, the wave speeds of different 
modes approach each other. This effect causes many mode jumpings as pointed 
out by Cotton and Sal wen [18] (in the case of a rotating pipe) and numerical 
calculations become extremely difficult. Also the mode jumping becomes more 
severe as n becomes large. Leibovich and Stewartson's analysis is valid only 
for the case of the most unstable mode and 1/2 q < «/n < 1/q with a/n suffi- 
ciently distant from the end points of the intervals. Finally, for a columnar 
vortex, they were able to obtain a sufficient condition which states that the 
flow is unstable if 


V 


dQ 

<Jr 


[ 


dQ dr 
HF cF 



(1.7) 


where Q is the local angular velocity and V is the velocity component in the 
azimuthal direction. Using a Taylor series expansion of the flow variables 
about a local radius r = r Q , Emanuel [28] showed that when viewed from a plane 
which contains all of the shear stress components, condition (1.7) is nothing 
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but a restatement of Rayleigh's criterion (centrifugal instability or inertial 
instability). In this particular plane then, the disturbances are helically 
symmetric. This idea which was first put forward by Pedley [14] has gained 
widespread attention in recent years. 

The results of [26] were extended by Stewartson [29] to include the 
effect of viscosity. Using asymptotic analysis, Stewartson was able to show 
that viscosity has a stabilizing influence which lowers the critical value of 
q, below which all disturbances are stable; and the most unstable mode occurs 
at a finite value of n. He found that when the Reynolds number Re >> 1, the 
sufficiency condition (1.7) for instability holds and the most unstable mode 

O /C 

has -n = 0(Re ' ). It must be emphasized that the results of [29] are valid 
only in the limit of high Re. 

Employing an asymptotic expansion, Stewartson and Capell [30] attempted 
to delineate the mode-jumping and separate the higher modes in the case of 
large n and a/n 1/q. They indicated (although not very convincingly) that 
u ->■ 0 as a/n + 1/q for all modes. The analysis of [30] breaks down 
completely for moderate values of the azimuthal wavenumber. In a follow on 
paper, Stewartson and Brown [31] showed the existence of near neutral center- 
modes (modes that are confined to the vicinity of r = 0) for moderate values 
of the azimuthal wavenumber, n. They reported that this instability exists 
only in distinct intervals of the swirl parameter, q. As n changes, so do 
the intervals. Since their analysis is valid for moderate values of n, it can 
and should be checked with numerical calculations. 

The conclusions of Leibovich and Stewartson [26] were confirmed by Duck 
[32] using a rather different asymptotic expansion. He showed a more complete 
description of higher unstable modes and his results indicate that the neutral 
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state is not attained but rather approached asymptotically. He obtained 
results for negative n as large as -15, which are in good agreement with his 
numer ical ly conducted cal cul at ions . 

It is hoped that this short literature survey has conveyed how difficult 
it is to obtain the stability characteristics of swirling flows -- even in the 
context of linear theory. Furthermore, it is the author's belief that the 
message conveyed by this survey is the fact that any future stability 
calculations should include the effects of viscosity. This was an important 
consideration in the formulation of the problem in this study and it will be 
shown later how fruitful the inclusion of viscous effects has been. 
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Chapter 2 


FORMULATION OF THE STABILITY PROBLEM 

A normal mode analysis is chosen for the stability calculations. In this 
approach, the basic flow state is perturbed by infinitesimal disturbances and 
the governing equations are then linearized with respect to these 
disturbances. The perturbations are decomposed into different Fourier modes 
and, because the equations are linear, each mode can be studied 

individually. It is further assumed that the flow is locally quasi-parallel 
and the mean velocity profile depends on the radial coordinate only. This 
last assumption is justified based on the fact that there is no external 

pressure gradient acting on the trailing line vortex. Viscous forces are 

responsible for any changes along the vortex axis and for high Reynolds number 
flow these changes are small over a wavelength. These assumptions are 

explained in more detail in the following sections. 

2.1 Mean Velocity Profile 

The similarity solutions for rotating porous pipe flows due to Donaldson 
and Sullivan [33] are selected as the mean velocity profiles. Their computed 
profiles, which are exact solutions to the three-dimensional Navier-Stokes 

equations of motion, show many of the flow features which are of interest in 
the study of unconfined trailing line vortices. For example, the solutions 
range from those which can be characterized as a single cell vortex to 
multiple cell vortices. In addition, experimental measurements have 
documented the existence of flows with many of these characteristics (see 
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Adams and Gilmore [34], Leuchter and Solignac [35], and Graham and Newman 
[36]). 

Cyl indr ical -pol ar coordinates (r, 9, z) are chosen as the coordinate 
system. Following Donaldson and Sullivan [33], the basic flow is taken to be 
of the form 


U = U(r) 

V = V(r) (2-1) 

W = z W(r) 

where, U, V and W are the radial, tangential and axial velocities, respec- 
tively. Note that the flow undergoes linear acceleration in the streamwise 
direction which is required from continuity considerations. 

The z-dependence of the axial velocity, W, indicated above is treated 
using a quasi-parallel flow approximation. The approximation can be justified 
on dimensional grounds using arguments similar to those employed by Donaldson 
and Sullivan in their similarity analysis. Assuming that the axial coordinate 
(z) has been made dimensionless by tube radius, R 0 , and using a Reynolds 
number. Re, based on tube radius, the axial velocity can be represented as 

W(r , z) - 4 C f'(n) (2* 2 ) 

r* . 

where n = C — ~ is a similarity variable for the radial coordinate, C 

R o 

represents a non-dimensional pressure gradient parameter, and r* is the 
dimensional radial coordinate. 

The maximum value for f'(ri) in Eq. (2.2) is unity. Consequently, if a 
reference location, z Q , is chosen as 


21 


then 


(2.3) 


z 0 = Re/4C, 

W(r, z Q ) = f'(n) • (2.4) 

At z q , the axial gradient of W(r, z) is given by 

I? V = R§ f ' ^ = * ( 2 - 5 ) 

which is very small since flows of interest have high Reynolds number and C is 
on the order of 0(1). The quasi-parallel flow assumption in the vicinity of 
z Q is therefore justified. For cases where z is small or where pressure 
gradient rather than viscous forces are responsible for the z-dependence, one 
can use a multiple-scale analysis to account for nonparallel effects (see for 
example [37]). However, in this study, only the case of a 'quasi parallel 1 
basic flow is presented. 


2.2 Governing Equations 

The governing equations of motion in nondimensional form for a steady, 
three-dimensional, constant density, Newtonian fluid (in cylindrical 
coordinates) are written as 


continuity 


1 1- (ru) +1-^ + .^ 
r dr v ; r 59 dz 


( 2 . 6 ) 


r-momentum 


5u 
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r 2~ 
[v u 


2 5v-i 

7aI J 


(2.7) 
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9-momentum 


1 dp 

at u 5r r 99 w dz r ' r ae 


+ ^ Tv^u v + ^ ® u_ | 

+ T?e Lv v ‘ "7 7? ae J 


( 2 . 8 ) 


and z-momentum 


aw . ~ aw . v aw , ~ aw _ ap , l „2~ 
af u ar r r ae w az" "azTCe vw 


(2.9) 


where 

2 2 2 

_2.a ,15 ,1 a ,a 

= 7^7 

Here, the nondimensional ization has been done with respect to the pipe radius, 
R 0 , centerline axial velocity, W 0 , time scale, R 0 /W 0 , dynamic pressure 

O 

magnitude p Q Vr, and Re is the Reynolds number based on R 0 and centerline 
velocity. The flow variables are assumed to consist of a mean part and an 
infinitesimally small perturbation, i.e. 


u = U + u 


v = V + v 


w = W + w 


( 2 . 10 ) 


p = n + p . 

Substituting Eq. (2.10) into the equations of motion (Eqs. (2. 6-2. 9)) and 
neglecting the quadratic terms, we obtain the linearized small disturbance 
equations 
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continuity 


r-momentum 


au + u + l5v.aw_ n 
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( 2 . 11 ) 
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i ^u 
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ae J 


( 2 . 12 ) 


0-momentum 
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(2.13) 


z-momentum 


aw 

af 


U SW + u^U^ + 
u ar u ar r ae 




_ ^E 
az 


+ i_ r $J!L + I aw 4 . 1 a 2 w . a 2 w n 

Re V 


(2.14) 


Since the coefficients of the perturbations in the above equations are 

functions of r only, the disturbance quantities are assumed to have a helical 
wave form: 


{u, v, w, p) = (i F(r) , G(r), H(r), P(r)} e 1 ( az+ne-wt) (2.15) 

Here, F, G, H and P are the complex disturbance eigenfunctions, a is the 
wavenumber in the axial direction, n is the wavenumber in the azimuthal 
direction, and w is the temporal frequency. In general, both a and ware 
complex variables. For a single valued solution, n must be an integer or 
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zero. 


When n is a 


When n equals zero the disturbances are axi symmetric, 
positive or negative integer, asymmetric perturbations exist which represent 
different spiralling directions of propagation, depending on the sign of n/u>. 

If spatial stability is considered, then w must be real and a = a + ia^. 
On the other hand, for a temporal solution, a is assumed real and w = + 

iw.. In either case, the sign of the imaginary part indicates decay or growth 
of the disturbance. 

Substituting Eq. (2.15) into the governing equations (Eqs. (2.11-2.14)), 
the linearized forms of the equations for perturbation eigenfunctions are 
obtained. They can be written in the form: 


continuity 


F 1 + £ + + ctH = 0 

r r 


(2.16) 


r-momentum 


i F" 
~ Re 


♦ 1 [« F- ♦ 1 iF-lr- 


He 


<n 2 +l 


a 2 )] 
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i 2n 
Re r‘ 


2V 


] G + P' = 0 


(2.17) 


0-momentum 



1 
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Re ' 2 


a 2 )] 


G + 




Re r 


2" + !*] F ♦ M. 

2 r r 
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(2.18) 
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z-momentum 


' Re + " Re ~ >^ H ’ + [ -ia) + ^r~ + iaW + <j> 

2 

+ R¥ ^~Z + a )] H + i ■gjpr F + iaP = 0 (2.19) 

where primes denote differentiation with respect to the radial coordinate. 

2.3 Boundary Conditions 

The boundary conditions for a variety of flow configurations have been 
given in detail in KMA. However, for the sake of completeness, they will be 
explained here. The boundary conditions at the outer wall are very simple 
since the perturbations must be equal to zero. In terms of the 
eigenfunctions, these conditions are 

F ( 1 ) = 0 

G(l) = 0 (2.20) 

H( 1 ) = 0 . 

Due to the singular nature of the coordinate system, and because all 
physical quantities must be continuous and bounded on the centerline, some 
non-trivial requirements exist as r ■+• 0. There are essentially two ways of 
treating this singularity. In one method, the entire pipe is treated as the 
computational domain with -1 <r <1 . Hence, condition (2.20) is enforced at r = 
± 1. However, in order to remove the centerline singularity, some parity 
conditions must be enforced. The nature of these parity conditions and the 
method of applying them are discussed by Metcalfe and Orszag [15], Leonard and 
Wray [38], and Orszag and Patera [39]. The second method treats the radius. 
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0<r<l, as the computational domain and hence the centerline is taken to be a 
boundary. Then, the need arises for a new set of relations to be enforced on 
the centerline. These compatiblity relations, given previously by Batchelor 
and Gill [40], and formalized recently by KMA are 


1 i m aft _ 

r+0 59 ‘ 


( 2 . 21 ) 


where V is the total velocity vector. This limit represents boundedness and 
smoothness conditions on the solutions along the centerline. The second 
approach has been employed in the present study. 

In expanding this condition (Eq. (2.21)), we need to consider only the 
perturbation part of the velocity since the mean flow is independent of the 
azimuthal direction. Representing the perturbation velocity field by v, we 
have 


dv - d fug +ve + w e ) 
59 50 e r 0 z 1 


( 2 . 22 ) 


or 
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-> dp 

1 im 5v _ 5u > r 5v > + 
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JT e 0 

d "0 ♦ 
d0“ = ' e r 


(2.24) 

(2.25) 

(2.26) 
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Substituting from Eq. (2.15), we deduce 


J+0 = nF - G) e p + (IF + inG) e 0 + inH e z = 0 . (2.27) 

In order for the equality to hold, each component of the resultant vector must 
be zero. Summarizing, in the limit along the centerline (r = 0), we have 

nF + G = 0 

F + nG = 0 (2.28) 


nH = 0 . 


The above conditions depend on 

the value of the azimuthal 

wavenumber, n, such 

that 



If n = 0 

F ( 0 ) = G(0) = 0 


and 

H ( 0 ) must be finite 

(2.29) 

If n * ± 1 

F(0) + G(0) = 0 


and 

H(0) = 0 

(2.30) 

If | n | > 1 

F(0) = G(0) = H(0) = 0 

(2.31) 


In the case when |n| = 1, two of the conditions become linearly dependent. In 
KMA, another relation was deduced by enforcing the continuity equation on the 
centerline, resulting in 


2F'(0) + nG'(O) = 0 . (2.32) 

However, after carefully examining the radial and azimuthal momentum equations 
as r +■ 0, it can be shown that condition (Eq. (2.32)) is identical to the 
conditions 
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F ' (0) = 0 


(2.33) 


or 

G ' (0 ) = 0 . (2.34) 

Therefore, in the present study, either Eq. (2.33) or (2.34) is used as the 
appropriate boundary condition for |n| =1. 
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Chapter 3 


NUMERICAL METHOD 
3.1 Introduction 

Over the past two decades spectral techniques have emerged as primary 
methods for solving hydrodynamic stability problems. The efficiency and high 
degree of accuracy obtainable by using these techniques has been demonstrated 
by others (see references [41, 15, 42, 43]). Very recently, KMA have applied 
a spectral collocation (pseudospectral ) method to a wide class of viscous, 
incompressible swirling flows. In their method, the flow variables were 
expanded in terms of a truncated Chebyshev series. The global eigenvalues of 
the discretized system were then obtained by a generalized complex QZ [44] 
routine. They have shown that the resulting algorithm is robust and easy to 
implement while being efficient. The advantages of their method over similar 
schemes are twofold. First, it is easily extendable to the case of 
compressible fields, since the coefficients of the flow variables are always 
evaluated in the physical space. Second, the use of the collocation 
formulation greatly simplifies the implementation of any new basic flow or 
coordinate transformation. The present study employs the same method, but 
with one major modification pertaining to pressure boundary conditions. 

The set of governing equations (2.16-2.19) are sixth order. Therefore, 
the six boundary conditions derived in Eqs. (2.20-2.34) are sufficient for a 
well posed boundary value problem. However, a straight forward 
discretization, as employed in KMA, would result in a system which requires 
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two artificial pressure boundary conditions. In KMA Neumann conditions were 
prescribed on the pressure at tie two boundaries. This pressure gradient 
condition is well known whenever the Navier-Stokes equations are solved in 
primitive variables, and has been discussed by Orszag and Israeli [45]. The 
prescribed values are obtained by taking the inner product of the vector 
momentum equation with the normal to the boundary. In cyl indr ical -pol ar 
coordinates, this turns out to be the r -momentum equation in its original form 
evaluated at the boundaries. KMA have shown that this formulation is viable 
and produces eigenvalues as accurate as any other method. However, in recent 
years it has become common to stagger the pressure cells or grid and eliminate 
the need for the two artificial pressure conditions. This is the approach 
taken here and is explained in detail in the next section. 


3.2 Chebyshev Spectral Collocation for Staggered Grid 

The unique properties of Chebyshev polynomials have been known for many 
years [46, 47]. Chebyshev polynomials distribute the error evenly, exhibit 
rapid convergence rates with increasing numbers of terms, and cluster the 
collocation points near the boundaries [48-50]. 

The kth Chebyshev polynomial, T ^ ( r: ) , is defined on the interval (-1, 1) 

by 


T k U) = cos[k cos' 1 ^] . 


(3.1) 


Because the spatial range in the present physical problem is (0, 1), a simple 
transformation is made from the physical variable, r, to the Chebyshev 
variable via 


l = 1 - 2r 


(3.2) 


3 L 


where 


- 1 < l < 1 . 


If 


l = COS <t> , 


it can be seen that relation (3.1) becomes 


T k (5) = cosk<t> . 


(3.3) 


Next a staggered grid is constructed such that the velocity components 

and the three momentum equations are evaluated at the cell faces, while 

pressure and the continuity equation are enforced at the mid cell points. A 

graphical representation of this grid is shown in Fig. 3.1. The velocities 

are evaluated at the collocation points which are the extrema of the last 

J 

retained Chebyshev polynomial (T^(^)) in the truncated series and are defined 
by 


Cj = cos JL , j = 0, 1, .... N (3.4) 

where the centerline and outer wall boundaries, correspond to j equals 0 and 
N, respectively. 

An interpolant polynomial is constructed in terms of the values of the 
flow variable at the collocation points by employing a truncated Chebyshev 
series. As an example, we present expressions for F(^) since extension to 
the other velocity components is straight forward. An Nth order polynomial 
for F(s) is written as 

N 

FU) = l g.(e) FU ) (3.5) 

j=0 J J 
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Fig. 3.1 A schematic of the staggered 

spectral collocation grid showing 
locations where F, G, H and P are' 
evaluated. 
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where the interpolant g.(S) is 

J 




(-D jtl (i - j?) t;( 5 ) 

7 sL 

C. 

J 


(? - 5j > 


(3.6) 


with 


C Q = C N = 2, Cj = 1, (1 < j < N-l ) 


It can easily be shown that 


(5 k> * 6 jk ■ 


(3.7) 


where 6^ is the usual Kroneker delta. Next, the first and second derivatives 
are determined explicitly by differentiating relationship (3.5). The 


derivatives, evaluated at £., are written as 

J 


dF 

dC 


N 

= 1 *jk F k 

j k=0 JK K 


(3.8) 


d 2 F 


d£ 


N 

= l B., F. 
j k=0 Jk k 


(3.9) 


j = 0, 1, ..., N 


where F^ is simply F(^), A and B are the elements of the derivative 
matrices and are given [50] by 


C. (-1) 
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k+j 


fljk ' S - V 


(j * k) 
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and 


(3.11) 


B., = A. A 
jk jm mk 


Now the pressure and the continuity equation are evaluated at the collocation 
points £j + 1 / 2 which are the roots t N ^) and are given 


^j+1/2 


= cos 


(2J+D 

2N 


j = 0, 1, 


N-l 


(3.12) 


Note that this set of points does not^ include the two boundary points. 
Next, the pressure is represented with an interpolating polynomial of 
N-l (rather than N) degree. Therefore, we write 


N-l 


P(e) = I h U) P(5. +1/2 ) . 
j=o J 


(3.13) 


where the interpolant h^U) is given as 


to 


(-l) j (sin <t i+1/2 T N U)) 

N ^ - C j+1/2 1 


j = 0, 


(3.14) 


N-l. 


Again, it can easily be shown that 


h j ( ^k+l/2 ) = 6 jk * 


(3.15) 


At this time, two sets of interpolating matrices are needed to 
interpolate from the staggered points to the grid points and vice versa. That 

is 


V 1/2 ■ k f 0 V F * ■ 


(3.16) 


j - 0, 1, ..*, N- 1 , 
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and 


P. 

J 


N-l 


l 

k=0 


M P 

jk k+1/2 * 


j = 0, 1 N , 


where the elements of the matrix M* are given as 


(3.17) 


★ 



<-n J+k+I 0 - 5 f tl/2 ) 1/2 

C k N (5 j+l/2 ' V 


(3.18) 


j = 0, 1, ..., N-l , 
k * o, 1, ..., N , 


and that of M as 


M.. = 
jk 


(-i) j+k (1 - c 2 




k+1/2- 


, 1/2 


C k+l/2~^ 


(3.19) 


j = 0, 1, ..., N , 
k = 0, 1, ..., N-l . 


In order for M* to be a square matrix, an extra row (j = N) with zero elements 
is added. Similarly in the case of M, a column (k = N) is added which 

contains null elements. 


In the case of the derivatives there are two different ways to 
interpolate. In one method the explicit functions g- and h- are first 

differentiated and then evaluated at the required points. Hence, the 
derivative matrices can be constructed explicitly. An equally valid method is 
to employ Eq. (3.16). That is, if 


dF 


N _* 
j+1/2 k=0 3 


(3.20) 


j = 0, 1, N-l , 
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then 


_* * 

A = M. A . 
jk jm mk 

Here we have used a mix of the two methods. That is, if 


dP 


N-l dh k (^,) 

j = Jo Pk+1/2 ’ 


j = 0, 1, ...» N , 


then we employ 
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where the elements of E are given as 
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A null column (k = N) is added to make E a square matrix. 

Finally, if the scaling factor for transformation between the 
and computational domain is given as 


(3.21) 


( 3 . 22 ) 


(3.23) 


(3.24) 


physical 
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C _ dc 
j ‘ 3r 


jj 0, 1, . . . , N , 


then the first derivative matrix, A, in the physical domain may be written 


as 


A., = S. A., 
jk J jk 


(3.25) 


and 


* * 

A,, = M. A , , 
jk jm mk 


(3.26) 


with similar relationships holding in the case of the mid cell points. 

Employing all of the relations developed above, the governing equations 
in discretized form are 
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0-momentum: 
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z-momentum 
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k=0 


The boundary conditions are: 


at 5 = - 1, 


F( -1 ) = G(-l ) = H(-l) = 0 for all n 


at £ = 1 


If n = 0 


F( 1 ) = G(l) = 0 
H 1 (1) = 0 . 


(3.31) 


(3.32a) 
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If n = + 1 

F( 1 ) ± 6(1) = 0 



H(l) = 0 

(3.32b) 


F'(l) = 0 or G * ( 1 ) = 0 . 


If | n | > 1 

F(l) = 6(1) = H(l) = 0 . 

(3.32c) 


3.3 Numerical Scheme 

The above equations are rearranged so that they can be represented in the 
generalized eigenvalue format as 

\ , (3.33) 

where D and L are coefficient matrices obtained from the appropriate 
discretized differential operator, and \ is an eigenvalue. For the temporal 
stability calculations, 

X = 0) , 
and 

$ = [F G H P] T . (3.34) 

Both D and L are square matrices with dimensions of 4N+3. In the case of 
temporal stability the last six rows of matrix D contain the boundary 
conditions. Since the boundary conditions do not contain the eigenvalue oj, 
the submatrix in D containing these conditions was made upper triangular 
through column operations. If all of the diagonal elements of this submatrix 
are non-zero, it indicates that the boundary conditions are independent of 
each other and the ranks of matrices D and L are reduced to 4N-3. 

It should be noted that the eigenvalue coefficient matrix, L, is 
singular. A procedure may be devised, using row and column operations, which 
reduces the rank of the coefficient matrix and removes the singularity (see. 
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e.g., Metcalfe and Orszag [15]). An alternative approach is to introduce a 
term pwP . in the continuity Eq. (3.27). This additional term makes the 
coefficient matrix for w non-singular and may be termed an artificial 
compressibility factor (see Malik and Poll [51]). Both methods were employed 
in KMA and it was found that computations using the method of artificial 
compressibility of Malik and Poll were at least one and one half times faster 
than the matrix operations of Metcalfe and Orszag [15], and both produced 
identical eigenvalues. Therefore, only the artificial compressibility 
approach is used in the present work. 

The artificial compressibility parameter, p, is assigned a very small 
value (on the order of 10' 18 ). This term generates large values for some of 
the eigenvalues of the matrix L"^D; however, experimentation with the value of 
p demonstrated that for values less than 10 ^ its effect on the desired 
(physical) eigenvalues is negligible. 

The general complex eigenvalue solver employed in the global calculations 
is the IMSL QZ routine called EIGZC. If the eigenfunctions are desired the 
eigenvalues are refined and the eigenfunctions are obtained subsequently using 
a local method. The local method which was employed here uses an inverse 
iterative technique (see Wilkinson [52]). The discretization of the governing 
equations for the local method is spectral but non-staggered. It was found 
that the eigenvalues obtained by the local procedure always converged to eight 
or nine significant digits within four iterations. Furthermore, it was found 
that both methods produced eigenvalues which agreed with up to six digit 
accuracy. 

All of the numerical solutions reported here were obtained on the CDC 
CYBER 860 and CYBER 205 machines at NASA Langley Research Center. 
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Chapter 4 


TEST CASES 
4.1 Introduction 

The convergence and accuracy of the non-staggered Chebyshev spectral 
collocation method has been tested extensively through numerous test cases in 
KMA. However, it is still desirable to obtain the convergence behavior of the 
staggered representation to contrast it with the non-staggered or tau 
formulation. The viability of the present method is examined here through two 
widely different flow configurations; namely Poiseuille flow in a pipe 
(confined flow) and for a Batchelor's type vortex (unconfined flow). In each 
case, numerous results have already been reported by others which greatly 
simpl ifies the comparison task. It must be emphasized that the objective here 
was not to duplicate all of the previous work, but rather to show briefly the 
efficiency and accuracy of the staggered collocation formulation. However, 
for one of the flows (Batchelor's vortex) where some quantitative differences 
with previous results were obtained, it will be necessary to explain the 
discrepancies. 


4.2 Temporal Stability of Poiseuille Flow in a Pipe 
The mean velocity for this problem is the well known profile given by 
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U = 0 


V = 0 

U = 1 - r 2 . 

The linear stability of Poiseuille flow in circular pipes subjected to 
azimuthally varying disturbances has been studied by Lessen, Sadler and Liu 
[53]. Using a shooting method, they found no instability for the n - 1 
perturbations. Two of their eigenvalues are tabulated in Table 4.1, along 
with the results of the present calculations. The results from the two 
methods agree very well. Next, the variation of the least stable mode with 
Reynolds number was obtained. A comparison of these results with those of 
Sal wen. Cotton and Grosch [54] is given in Table 4.2. The agreement is 

excel lent. 

The convergence behavior of the present method as compared to the 
Chebyshev tau formulation of Metcalfe and Orszag [15] and the non-staggered 
method of KMA is presented in Table 4.3. It is clear that Chebyshev 

collocation (in either form) has a much better convergence rate than that of 
the tau method. However, there is hardly any difference between the 

convergence rate of the two collocation formulations. This fact suggests 
that, contrary to prevailing opinion, when the artificial pressure boundary 
conditions are employed correctly they do not affect the accuracy or the 
convergence rate of the spectral methods. 

Although not reported here, the two collocation methods have been 
compared for a variety of other flows. In each case, both methods produced 
identical results up to seven or eight significant digits. 
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Table 4.1 Comparison of the First Eigenvalue Predicted by 
Lessen et al . [53] with those Using Chebyshev 
Collocation 





Lessen et al . [53] 

Present 

Method 

n 

a 

R e 

C r C.j 

C r 

Ci 

i 
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200 

0.645 -0.129 

0.64526 

-0.129205 

i 

i 

2200 

-0.067 

0.39797 

-0.067709 

i 

i 

2200 

— — 

0.89663 

-0.048114 
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Table 4.3 The Convergence Behavior of the Two Least Stable 

Modes for Poiseuille Flow in a Pipe. N is the number 
of polynomials used to resolve each flow variable. 
Here, Re = 10, a = 1 , and n = 1 . 


Metcalfe and Orszag [15], Chebyshev Tau Formulation 


0.501654048- i 1 . 392398524 
0. 491070208- i 1 . 393534481 
0. 491063984- i 1.393490921 
0. 491064084- i 1 . 393490894 

Present Calculations, Staggered 

0. 491065536- i 1.3934907 13 
0. 491064085- i 1.393490894 
0. 491064084- i 1.393490894 
0. 491064084- i 1.393490894 

KMA, 

0. 491067022- i 1.393495866 
0. 491064084- i 1 . 393490896 
0. 491064084- i 1.393490894 
0. 491064084- i 1 . 393490894 


0J2 

0. 780474345- i 2. 745100739 
0. 763978453- i 2. 807040555 
0. 762026865- i 2. 807291233 
0.76202421 5- i2. 807286430 

Chebyshev Collocation 

0. 762034238- i 2. 807291 149 
0. 762024225- i 2. 807286425 
0. 762024223- i 2. 807286422 
0. 762024223- i 2. 807286422 


0.762035973- i 2. 807301022 
0. 762024226- i 2. 807286427 
0. 762024223- i 2.807286422 
0.762024223-i2. 807286421 


Non-Staggered Chebyshev Collocation 


For Poiseuille flow, 16 Chebyshev polynomials were required for eight figure 
accuracy. The number of polynomials, N, depends on the particular mean flow, 
and between 32 and 70 polynomials were used in the present study to insure 
convergence. 

4.3 Temporal Stability of a Trailing Line Vortex 

The approximate asymptotic form of the solution for the mean velocity 
profile of a Batchelor's vortex [24] is 


U = 0 

V = -p — [1 - exp(-r^)] (4.2) 

W p 

W = ~ + exp(-r ) . 
w 2 

where the azimuthal and axial velocities, have been normalized with respect to 

W 2 » which is the difference between the centerline velocity and the far field 

velocity. The quantity, q, in Eq. (4.2) has been called the swirl parameter 

by Lessen et al . [24] and is related to the ratio of the maximum swirl 

velocity to the maximum axial velocity excess (or defect) • Throughout the 

W 1 

present study, it is assumed that rr- - 1 and therefore 

2 

W = 1 + exp(-r^) . (4.3) 

The two components of the mean velocity profile are shown in figure 4.1. It 
has been pointed out by Lessen et al . [24] that the uniform part of 


*If q is the ratio of maximum swirl velocity to the maximum axial velocity, 
then q = 0.639 q in this case. 
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1.0 


a) 



Fig. 4.1 Velocity components for a Batchelor 
vortex. 
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the flow in Eq. (4.2) only affects the frequency and not the growth rate of 
the disturbance. 

Since r varies between zero and infinity, it is necessary to map the 
coordinate system into a Chebyshev (computational) domain. An algebraic 
transformation given by Malik, Zang and Hussaini [55] is employed in the form 


r 


1 + l 


(4.4) 


where a and b are constants. For any arbitrary value of a, b is uniquely 
determined by 


It should be noted here that parameter, a, is very significant and that the 
above transformation concentrates at least half of the collocation points 
between r = 0 and r = a. Throughout this study, the values of a and r max were 
fixed at 3 and 100, respectively. The convergence behavior of the transforma- 
tion (Eq. (4.4)) is fully explained in Grosch and Orszag [56] and KMA. Also 
the tests conducted for the selection of a and r max have been discussed in 
detail in KMA and will not be repeated here. 


Figure 4.2 shows the variation of the growth rate of the first three 
unstable modes with the axial wavenumber, a, for n = -2, and q = 0.8, at a 
Reynolds number of 10,000. The curves are identical to the ones resulting 
from the inviscid calculations of Duck and Foster [19]. Since the present 
viscous calculations produce identical eigenvalues to those of Duck and 
Foster, it is logical to assume that the instability is an inviscid one. This 
would imply that the number of unstable modes depends directly on the number 
of collocation points as noted by Duck and Foster. However, that conclusion 
does not appear to be consistent with experiments, and further results in the 
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present study will shed some light on this problem. As pointed out by [19], 
finer resolution is needed to resolve higher modes. Note that the first mode 
corresponds to that found by Lessen et al . [24]. The variation of the real 
part of the frequency (oo r ) vs a is shown in figure 4.3. The curve 

corresponds to mode 1 in the previous figure. Since the other two modes have 
real parts which are nearly identical with the one shown, they are not plotted 
for clarity. The variation of growth rate vs axial wavenumber for the case of 
n = -1 is shown in figure 4.4. The calculations were made at Re = 10,000, q 
0.4 and only results for the first three modes are plotted. Again, the 

primary mode corresponds to the inviscid mode obtained by [24]. The cor- 
responding to this mode is plotted in figure 4.5. Next, the variation of the 
first two modes with Reynolds number for the case of n = - 2 at two different 
values of a are shown in figures 4.6 and 4.7. It is apparent from these plots 
that they are indeed inviscid modes for Reynolds numbers above 10,000, since 
they have almost attained their maximum growth rates and are virtually 
independent of Re. In the case of n = +2, no instability has been observed. 
Although the results shown in figures 4.8-4.10 confirm this fact, figure 4.10 
indicates that as Reynolds number increases, the flow approaches a marginally 
stable condition for this azimuthal wavenumber. 

Even though the Batchelor vortex was chosen primarily as a test case to 
evaluate the staggered collocation method, new insights and some fundamental 
results have been produced. The fundamental results are related to the 
identification of new viscous instability modes and were produced during the 
testing. Before presenting the discussion of the new results, as mentioned in 
Sec. 4.1, it is important to explain some discrepancies which exist between 
this study and earlier work. 
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Fig. 4.3 Variation of the real part of 

frequency of n = -2 disturbances 
with wavenumber^for a trail ing ’line 
vortex. Here, 1 - 0*8, Re = 1(F, 
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Fiq. 4.4 Variation of the qrowth rate of 
asymmetric (n =-l) disturbances 
with wavenumber^for a trailing 
vortex. Here, q 3 0.4, Re = 10 








Fig. 4.7 Variation of the real part of 

frequency of n = -2 disturbances 
with Reynolds number for a trailing 
line vortex. Here, ci = 0.8. 
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9 Variation of the real part of 
frequency of n = +2 disturbances 
with wavenumber for a trailing line 
vortex. Here, q = 0.8, Re = 10 4 . 
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Recently, results were reported in KMA which indicated that although 
qualitative agreement exists between the collocation method and the viscous 
calculations of Lessen and Pa ill et [27], there are quantitative disagreements. 
The disagreement is shown in figure 4.11 (taken from KMA) which presents the 
variation of u^. , with the swirl parameter q. Of interest is the existence 
of a second mode of instability which was not obtained by [27]. The staggered 
collocation method employed here has confirmed KMA's earlier findings and 
their two curves shown in the figure were reproduced identically. The above 
result indicates that these inviscid higher modes exist at much lower Reynolds 
numbers than had been expected. A further test was conducted in which the 
critical values, given by Lessen and Pail let [27], were used as input for the 
present method and was calculated. The results are tabulated in Table 
4.4. It is clear from the magnitude of w/s that the values taken from [27] 
are not very close to the critical point. In order to show the accuracy 

obtainable by the collocation method, the critical value for the axisymmetric 
case is presented in Table 4.5. After careful examination, it was determined 
that the discrepancies are only with the viscous results of [27]. The 
inviscid calculations of Lessen et al . [24], which employed similar methods to 
those used in the later viscous analysis, are consistent with the collocation 
results. The discrepancy can be explained by a detailed examination of the 
approach taken by Lessen and Pail let. 

Starting the solution at r = 0 with a Frobenius series and at r = r max 
with an asymptotic solution. Lessen and Paillet [27] integrated from both 
limits using a Taylor series expansion and matching the solutions at some 
intermediate radius. They mentioned that the results were strongly dependent 
on the radius of integration, r max as well as the intermediate matching 
condition. The dependency of the eigenvalue on r max has been demonstrated by 
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Results of Lessen and Paillet 
KMA study, first mode 
KMA study, second mode 


Fiq. 4.11 Variation of the cjrowth rate with 
9 swirl parameter, q, for a trailing 

line vortex. Here, n - -2, 
n = 1.34 and Re = 141.4 • 




Table 4.4 The Critical Values of w. Obtained by the Present 
Method for Different Azimuthal Wavenumbers. The 

Critical Values for Re , a , q were Taken from Lessen 
and Pail let [27]. ~ C 


n 

a c 

<c 

Re c 

“i 

-i 

0.42 

0.45 

13.9 

-0.002547 

-2 

0.91 

0.70 

27.9 

-0.029787 

-3 

1.62 

0.95 

48.2 

-0.05455 


Table 4.5 Critical Values of Different Parameters 
for the Axisymmetric Case n = 0. 


a c 

Re c 

% 

00 

0.468 

322.35 

1.08 

0.441679- i 4. 97x1 0 -7 
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further testing. Table 4.6 shows the variation of the maximum growth rate for 
different azimuthal wavenumber, n, as compared with the inviscid calculations 
of [19, 24]. Remember that the values of Re apply only to the present viscous 
calculations. The results obtained by the present method indicate very good 
agreement with those of [19, 24]. In the case of n = - 1, as the Reynolds 
number increases, the convergence of the eigenvalue toward the value obtained 
by the finite difference global method of Duck and Foster [19] is excellent, 
while it is acceptable with respect to that of [24]. Duck and Foster [19] 
have indicated that for n = -1, they had to go to larger values of r max such 
as 10 or 12 in order to obtain convergent results. The present result was 
obtained with r max = 100, which is much larger than their r max , and 
substantiates their findings. Finally, the dependency of the eigenvalues on 
r max clearly demonstrated in the work of Paillet [57] and was recently 
confirmed by Lessen (private communication). 

It is well known that in the previous stability calculations of the 
Batchelor vortex, axisymmetric instability was not observed. However, while 
conducting some routine calculations in the present study, two new modes of 
instability were discovered. One instability mode was an axisymmetric, n = 0, 
mode and the other was an asymmetric, n = 1, mode. Both are viscous modes. 
Figure 4.12 shows the variation of growth rate, ^ with a for n = 0, q = 1 
and Re = 10,000. The figure indicates that this mode has a growth rate which 
is orders of magnitude smaller than the inviscid modes. The maximum value 
of u).j occurs roughly at a = 0.3 which suggests that the disturbance's 
wavelength is on the order of the core radius of the vortex. The variation of 
the real part of the frequency with a is plotted in figure 4.13. The nearly 
straight line indicates that the axial phase speed of the wave, C r = uy/ a, is 
constant and equal to unity over a wide range of a. Variation of w vs 
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Table 4.6 Comparison of Maximum Growth Rate for 
Different Azimuthal Wavenumber. 


Re 

n 

a 

max 

/V 

%iax 

Inviscid 
Lessen 
et al . 
[24] 

w i 

Inviscid 
Duck and 
Foster 
[19] 

w i 

Present 

Method 

“i 

20000 

-1 

0.3 

0.32 

0.1470 

0.14955 

0.14904 

100000 

-1 

0.3 

0.32 

0.1470 

0.14955 

0.14935 

20000 

-2 

1.2 

0.70 

0.3138 

0.31392 

0.31274 

50000 

-2 

1.2 

0.70 

0.3138 

0.31392 

0.31332 

100000 

-2 

1.2 

0.70 

0.3138 

0.31392 

0.31352 

100000 

-3 

1.7 

0.79 

0.3544 

0.35462 

0.35404 

100000 

-4 

2.15 

0.82 

0.3777 

0.37754 

0.37687 
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Reynolds number are presented in figures 4.14 and 4.15, for two distinct 
wavenumbers. As the Reynolds number increases, the disturbance attains its 
maximum growth rate very rapidly and then drops gradually. The maximum of w. 
shifts to higher Re as a increases and there is a significant drop in its 
magnitude. Figure 4.15 shows that the phase speed is constant over a wide 
range of Reynolds number and has a magnitude which is slightly less than the 
uniform outer flow. The variation of the wavenumber of the maximum growth 
rate with Reynolds number is shown in figure 4.16. The curve is consistent 
with the behavior of a viscous mode and shows the damping effect of viscosity 
on the wavenumber at moderate Reynolds numbers. Finally, the critical values 
are compiled for the different control parameters for this axisymmetric mode 
in Table 4.5. 

The axisymmetric mode just described has been observed in many 
experiments and photographs of many contrails of commercial aircraft. 
Although the phenomena is termed "core bulging" or "core bursting" in the open 
literature, the present author believes that is a misnomer since the phenomena 
is traveling outside the core and the core remains intact. This fact is 
demonstrated in the photographs shown in figures 4.17-4.19. These photos 
which are contrails of a large commercial aircraft were reproduced from the 
work of Bisgood [58]. They were taken from the ground with the vortices at 
heights between 10,000 and 12,000 meters. The separation distance between the 
vortices was between 25 and 40 meters. Figure 4.17 shows clearly that, while 
axisymmetric instability waves are present, the core still remains intact. 
The next two photos (figures 4.18 and 4.19) show the amplification and growth 
of these waves. It must be emphasized here that these photographs are not 
sequential and might not even be the same contrails. However, of importance 
here is the presence of the Crow [59] instability in every photo. It seems 
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a = 1.0 


!_[ 


4.15 


O.i 



10 4 


Re 


Variation of the real part of 
frequency of axi symmetric (n = 0) 
disturbances with Reynolds number 
for a trailing line vortex. 

Here, q = 1.0 • 
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that the two forms of instability (core bulging and Crow) do not interact 
significantly. This fact was established recently by Sarpkaya and Daly [60] 
who found that any combination of the instabilities might occur. 

Figure 4.20 is another photograph which has been reproduced from the work 
of Bisgood [58]. Sarpkaya's finding can be supported by this picture which 
shows Crow's classical "sinuosities" mode without any sign of the axisymmetric 
wave form. Until further research is carried out, further comments on the 
nature of the interaction between these two distinct forms of instabilities 
are speculative and therefore unjustified. 

The core radius where the azimuthal mean velocity reaches its maximum is 
approximately at r = 1.1. In order to substantiate the fact that the 

axisymmetric wave is occurring just outside of the core of the vortex, one has 
to look at the perturbation eigenfunctions. Figures 4.21a-4.21h represent the 
eigenfunction amplitudes and phase angles for all of the perturbation 

components near their critical values. The eigenfunction amplitudes have been 
normalized by the maximum amplitude of the azimuthal eigenfunction throughout 
this study. Figures 4.21a and 4.21c are significant and show clearly that 

most of the perturbation energy in the radial and azimuthal directions is 
spread over a wide distance and is mostly outside of the core. However, there 
is hardly anything happening beyond r = 10. The phase plot for G in figure 
4.21d suggests the existence of a critical layer just outside of the core at r 
1.5. A similar critical layer exists in the case of the axial perturbation 

(figure 4.21f); however, it is positioned inside of the core radius at a 

distance of about r = 0.7. This is understandable since the mean axial 
velocity has an inflectional point close to r = 0.75. 
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Fig. 4.20 Crow's [59] classical 
"sinuosities" mode . 



Amplitude 
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Fig. 4.21 Spatial distribution of amplitude and phase of various 
eigenfunction components of an axisymmetric (n = 0) 
perturbation for a trailing line vortex, a) and b) radial 
velocity, c) and d) azimuthal velocity, e) and f) axial 
velocity, g) and h) pressure. Here, a = 0.4. a = l o 
Re = 400. * M 
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Fig. 4.21c) Continued. 


78 


Phase 



79 




r 

Fig. 4.21e) Continued. 
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Amplitude 



Fig. 4.21g) Continued. 
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The variation of the growth rate, a,., with the axial wavenumber, «, for 
the asymmetric mode n = 1, is plotted in figure 4.22 with q = 0.4 and Re = 
10,000. Of the two modes shown, the second mode is stable. The first mode 
has a growth rate comparable to the axisymmetric wave but is unstable over a 
smaller range of a. However, like the n = 0 case, it is basically a long 
wave instability. The variation of maximum growth rate of mode 1 with a is 
shown in figure 4.23. The maximum occurs at a . 0.34. The variation of u> r 
with « for the first mode is presented in figure 4.24. The nearly linear 
variation indicates that the axial phase speed of the wave is slightly greater 
than unity and, like the axisymmetric case, is traveling with the uniform 
outer mean flow, although it is a bit closer to the core of the vortex. This 
form of instability has been observed experimentally both in wind tunnel 
measurements and contrail photos. The measurements of Singh and Uberoi [61] 
and Strange and Harvey [62] support the existence of an asymmetric disturbance 
which has an axial wavelength on the order of the core radius. Furthermore, 
the photograph of condensation trails behind a wide-bodied jet shown in [62] 
(their Fig. l) clearly displays the existence of helical disturbances which 
co-exist along with the Crow instability. The viscous nature of the n = 1 
mode is shown in figures 4.25 and 4.26. The variation of growth rate w . 
with axial wavenumber, a, shows similar behavior to that of the n = 0 
disturbance. However, in this case, the instability occurs at much lower 
Reynolds numbers. Also the maximum growth rate which occurs at Re « 60 is an 
order of magnitude greater than the . imax for the n - 0 case. But the most 
interesting phenomenon is the behavior of these two instabilities at high 
Reynolds numbers. While the asymmetric mode is the stronger of the two 
instabilities at low Reynolds numbers, as Re increases above 10 4 , the 
axisymmetric disturbance becomes the mode with the higher growth rate. This 
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Fiq. 4.26 Influence of Reynolds number on 
the real part of frequency of 
asymmetric (n = +1) disturbances 
for a trailing l|ne vortex. 

Here, a = 0.34, q = 0.4 . 
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may be why, in most of the contrails behind large commercial jets (which occur 
at high Reynolds numbers), the axisymmetric mode is the dominant and 

persistent form, while in wind tunnel tests, which are conducted at much lower 
Reynolds numbers, both modes are present equally. The amplitudes and the 

phases of the disturbance eigenfunctions are presented in figures 4.27a-4.27h. 
Remember that disturbances having |n| = 1 are the only ones with non-zero 

radial and azimuthal velocity on the centerline. This fact is clearly 
demonstrated in figures 4.27a and 4.27c. Although the peaks occur on the 

centerline, a considerable amount of the perturbation's energy is just outside 
of the core. This feature was explained earlier when discussing the axial 
phase speed of the wave. Figure 4.27e shows the variation of the perturbation 
amplitude in the axial direction with radial distance. The peak is positioned 
just inside of the core with a significant part of the perturbation imme- 
diately adjacent to the outer part of the core. The overall observation of 
the eigenfunctions indicates that this instability is much tighter in radial 
extent than that of the n = 0, disturbance. A significant point to be made 
here is the fact that, unlike inviscid modes, no higher unstable modes were 
found for either the axisymmetric or asymmetric form of instability. This was 
tested extensively by going to successively higher values of N (number of 
Chebyshev polynomials) than the number needed to obtain convergence. This 
point has to be explored further. 

Much work needs to be done with regard to documenting these two new 
viscous modes of instability. Of particular importance is the variation of 
the perturbation quantities with respect to the swirl parameter q. However, 
neither space nor time allows for such an extensive study here (Recall that 
this was only to be a test case for calibrating the numerics.). A more 
detailed study should follow in the near future and will be reported 
el sewhere. 
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Amplitude 



Fiq. 4.27 Spatial distribution of amplitude and phase of various 
eigenfunction components of an asymmetric (n = +1J 
perturbation for a trailing line vortex, a) and b) radial 
velocity, c) and d) azimuthal velocity, e) and f) axial 
velocity, g) and h) pressure. Here, a = 0.34, q = 0.4, 
and Re = 50. 
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Fig. 4.27c) Continued. 
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Fig. 4.27d) Continued. 
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Fig. 4.27e) Continued. 
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Fig. 4.27g) Continued. 
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Finally, before concluding this chapter, two important facts related to 
the discovery of the new viscous modes are: 

1. Batchelor's vortex represents a good approximation to 
trailing line vortices and can be employed for further 
stability studies. 

2. Although the stability calculations were conducted for a 
laminar vortex, the excellent agreement obtained with 
experimental measurements and contrail photos, which 
involve mostly turbulent vortices, indicate that these 
long wave viscous instability modes, once activated, do 
not interact significantly with the surrounding turbulent 
field and are not influenced strongly by the magnitude of 
the scaling parameters. 
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Chapter 5 


RESULTS AND DISCUSSION 

In this chapter, the results obtained using the three-dimensional 

profiles of Donaldson and Sullivan [33] are presented. Due to the wealth and 

variety of solutions possible from their compilation as well as constraints on 

computer resources and time, representative velocity profiles had to be 
chosen. In each case, care was taken to insure that the selected profiles 

approximated trailing line vortices as closely as possible. In section one, 
the stability of a single cell vortex is studied. The second section contains 
the results for a comprehensive study of a two cell vortex. The third section 
is devoted to the case of a three cell vortex which is studied less 
extensively than that of a two cell vortex. Except in the case of the three 
cell vortex where the effect of high values of azimuthal wavenumber, n, on the 

stability characteristics were obtained, the calculations are limited to n = 
0, +1, and +2. 

Before proceeding further, two points need to be emphasized. First, 
unless specified otherwise, the axial gradient term, <p, as defined in Eq. 
(2.5), is included in the governing equations even though all of the 
calculations are conducted at high Reynolds numbers. It will be shown that 
although the effect of this term on the modes with negative n is negligible, 
it causes drastic changes in the case of positive azimuthal wavenumbers. 
Secondly, rather than using q (as defined in Eq. (4.2)) the alternate 
parameter, q, which is defined as the ratio of the maximum swirl velocity to 
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the centerline axial velocity, has been used to delineate the stability 
characteristics. 

As mentioned in the introduction, in the case of swirling pipe flows, the 
existence of center- as well as wall-modes of instability has been established 
by previous research efforts. In addition, previous investigations have 
identified ring modes as those which occur at some intermediate radius. This 
investigation has shown that compound instability modes which will be called 
co-instabilities occur for many types of swirling flows. The results of the 
present work show, for the first time, cases where center instability modes 
coexist with wall modes and/or with ring modes, as well as other 
combinations. These different types of compound instability modes will be 
described as they occur. It is obvious here that mode identification has 
become a confusing and complex matter. Originally, the word "mode" was used 
to separate disturbances having different azimuthal wavenumbers. The 
discovery of higher modes associated with each negative value of n and now the 
presence of center-, ring- and wall -modes of instability has created an 
untenable situation. Thus, to simplify the text (and stop repeating "mode"), 
from here on center modes will be indicated by C, ring modes by R and wall 
modes by W. Therefore a combination of a center mode and a ring mode is 

identified as C R with similar identification applying to any other 
combinations. The word mode will be used only in its original context, that 
is to represent a specific value of the azimuthal wavenumber and the higher 
modes associated with that particular value. 
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5.1 Stability of a Single Cell Vortex at Re = 4000 

The three components of the velocity profile associated with this vortex 
are shown in figure 5.1. The radial variation of the azimuthal velocity V, 
presented in figure 5.1b, is identical to that of an unconfined trailing line 
vortex. The variation of the growth rate with axial wavenumber for an 
axisymmetric disturbance (n = 0) is shown in figure 5.2 for q = 0.19 and 
1.0. No sign of instability was observed. Furthermore, as a increases the 
flow becomes more stable. The variation of co p with a is presented in figure 
5.3. It seems that for q = 1.0, the axisymmetric mode is almost stationary. 
This is in accordance with the criterion (Eq. (1.6)) obtained by Ito et al . 
[21]. The effect of Reynolds number on the axisymmetric mode at q = 1.0 is 
presented in figures 5.4 and 5.5 for a = 1.0 and 0.3, respectively. Although 

the flow becomes less stable as Re increases, an unstable state is never 
reached. 

The variation of the growth rate with a for the asymmetric mode (n = 1) 
is shown in figure 5.6 for q = 0.8 and 0.4. The flow is clearly stable and 
becomes more stable as a increases. The real part of frequency versus axial 
wavenumber is plotted in figure 5.7. The calculations for the case n = 2 
produced similar results, showing mainly that this profile is highly stable 
with respect to the modes having n > 0; at least for the range of parameters 
studied here. The above conclusion has been reached in the case of inviscid 
disturbances for Long's vortex by Foster and Duck [23]. 

Figure 5.8 shows the variation of the growth rate, « i# with a associated 
with the asymmetric mode n = -1, for q = 0.4. Similar to the case of a 

trailing line vortex, the instability is a very explosive one with higher 
modes present. Here, only the two most unstable modes are plotted for 
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Fig. 5.2 Variation of the growth rate of 

axi symmetric (n = 0) disturbances 
with wavenumber for a single cell 
vortex. Here, Re = 4000, * 0 • 





5.3 Variation of the real part of 

frequency of axi symmetric (n = 0) 
disturbances with wavenumber for a 
single cell vortex. Here, 

Re = 4000, * 0 . 
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Fig. 5.5 

Influence of Reynolds number on the 
real part of frequency of ^ 

ax i symmetric (n - UJ a 
?„r a single cell vortex. Here, 

q = 1.0, * * 0 • 



Fig. 5.6 Variation of the growth rate of 
asymmetric {n = +1) disturbances 
with wavenumber for a single cell 
vortex. Here, Re = 4000, u, * 0 . 



5.7 Variation of the real part of 

frequency of asymmetric (n - +1J 
disturbances with wavenumber for a 
single cell vortex. Here, 

Re = 4000 , il< * 0 • 







clarity. However, there are at least two major differences between this 
instability and the one associated with Batchelor's vortex. First, this 
instability exists at higher axial wavenumbers, a. Second, both stability 

modes show evidence of co-instabilities. That is, lobes are apparent on 
either side of a = 0.7. Varying a does not show co-instabilities directly, 
and in fact, co-instabilities could exist when growth rate curves show no 
evidence of unusual behavior. It is important to note that change near a = 
0.7, suggests a shift in the dominant type of instability. Furthermore, it is 
noted that the stability curves in this figure are similar to the trailing 
line vortex curves shown in figure 4.4, for 0 < a < 0.7. The co-instabil- 
ities were checked by looking at the eigenfunctions of the disturbance at 
several values of a. That study of eigenfunctions revealed that this indeed 
was the case, which pointed to a disturbance having a very complex structure 
with several peaks in its ampl itude. 

The possibility of coexistence of several different modes was first 
mentioned in the work of Leibovich and Stewartson [26]. Subsequently, it was 
shown in the inviscid asymptotic analysis of Duck [32] that, for large 
negative n, there are two critical layers associated with the Batchelor's 
profile. He showed that the peaks in the disturbance eigenfunction are 
situated near these critical points (incidentally his eigenfunction curves 
show this multiple peak phenomena, but only for the higher modes, with no such 
effect for the primary mode). In the case of small negative n. Duck has 
stated that there is only one critical layer in the flow field. Depending 
upon the value of a, multiple peaks in the eigenfunctions for the single and 
multiple cell vortices do exist in the present calculations. We will discuss 
this later in the context of two cell vortices. There we will show that the 
peaks in the eigenfunctions could be used to distinguish the type or types of 
modes present (namely center, ring and wall modes). 
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We will also discuss this co-instability feature in detail for the more 
interesting case of a two cell vortex. Setting 4 ; = 0, the calculations were 
repeated and the results are plotted in figure 5.9. A comparison of figures 
5.8 and 5.9 indicates that the effect of the axial gradient term, <|>, is 
negligible. Actually, this turned out to be true (in the case of single cell 
as well as multiple cell vortices) whenever n < 0. The changes in oj r with 
a, for the two modes shown in figure 5.8, are plotted in figure 5.10. Notice 
the kink in the higher mode curve which occurs roughtly at a = 0.7. The 
inviscid nature of the asymmetric mode, n = -1, is revealed in figure 
5.11 for q = 0.4 and a = 1.0. The behavior is very similar to the trailing 
line vortex case presented in figure 4.6. As the Reynolds number increases, 
the growth rates, w. , increase very rapidly and asymptote to constant 
values. In fact, in the case of the primary mode at Re = 4000, the inviscid 
value of the growth rate is almost achieved. The variation of the real part 
of frequency, w,, with Re is shown in figure 5.12. Variation in the real 
part of u over a wide range of Reynolds numbers is small. The effect of 
higher swirl on the growth rate of the asymmetric mode (n = -1) is presented 
in figure 5.13. At q = 0.8, the lobe in the higher mode has disappeared while 
primary mode "splitting" has become more acute. There is a substantial 
increase in the growth rate of the first mode, particularly at higher values 
of axial wavenumber, a. The variation of the real part of the frequency 
with a for this disturbance is shown in figure 5.14 for q = 0.8. Notice the 
shift in the w, of the first mode (which results in the kink) even though 
the two distinct sections of the curve have similar slopes. 

The behavior of the growth rate with a for disturbances having n = -2 is 
shown in figure 5.15 for q = 0.8. Here, the peaks associated with the lower 
values of the axial wavenumber are the dominant ones. Overall the growth rate 
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5 11 Influence of Reynolds number on 
the growth rates of asymmetric 
(n = -1) disturbances for a 
single vortex. Here, a = 1.0, 

= 0.4 and il» * 0 . 



Fig. 5.12 Influence of Reynolds number on 
the real part of frequency of 
asymmetric (n =-l) disturbances 
for a single cell vortex. 

Here, a = 1.0, q = 0.4, and ij* * 0 



5.13 Variation of the growth rate of 
asymmetric (n = -1) disturbances 
with wavenumber for a single cell 
vortex. Here, q = 0.8, Re = 4000 
and iji ^ 0 • 



Ficj. 5.14 Variation of the real part of 

frequency of asymmetric (n = -1) 
disturbances with wavenumber for a 
single cell vortex. Here, 
q = 0.8, Re = 4000, and p * 0 . 
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it i on of the growth rate of 
-2 disturbances with 
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indicates a twofold increase over the values of for Batchelor's vortex 

(figure 4.2). Also, the instability presists over a much larger range of the 
axial wavenumber, a. The problem of mode identification, as pointed out by 
[26], is demonstrated in figure 5.16. Note how closely the two curves follow 
each other. In fact, as the neutral curve is approached, the real part of 
both frequencies become indistinguishable when they are displayed 
graphically. However, contrary to the remarks made by Leibovich and 
Stewartson [26], at no time did we encounter difficulties with the convergence 
of the first few eigenvalues. Rather, the difficulty is caused by the 
successively smaller increments of a needed to follow the same eigenmode. 
This in turn would make global calculations very expensive rather than 
impossible. 

The variation of the growth rate with Re for the n = -2 mode at q = 0.8 
and a = 2 is presented in figure 5.17. Clearly the instability is an 
inviscid one. The real part of the frequency, shown in figure 5.18, indicates 
little variation with Reynolds number, especially at high Re. 

5.2 Stability of a Two Cell Vortex at Re = 10000 

The three components of the velocity profile are shown in figure 5.19. 
Again, note the similarity between the azimuthal velocity profile to that of a 
trailing line vortex. Also there is an inflectional point in the axial 
component of velocity between r = 0.6 and r = 0.8. The variation of growth 
rate with axial wavenumber for the axisymmetric mode (n = 0) is shown in 
figure 5.20, for q = 0.5 and 1.0. Obviously, increasing the swirl has a 

destabilizing influence. However, the surprising effect here is the enormity 
of the growth rates obtainable with this mode. In fact, the maximum growth 
rate is orders of magnitude larger than the viscous mode of the Batchelor's 
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5.16 Variation of the real part of 

frequency of n = -2 disturbances 
with wavenumber for a single cell 
vortex. Here, q = 0.8, Re = 400C 
and di / 0 . 



Fig. 5.17 Influence of Reynolds number on 
the growth rate of n = -2 
disturbances for a single cell 
vortex. Here, a = 2.0, q = 0.8 
and ip * 0. 




Ux 10 





Fig. 5.19 Velocity components for a two cell 
vortex. 
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vortex. No axisymmetric mode of instability has been reported previously 
(regardless of whether the vortex was single cell or multiple cell) and the 
explosive growth rates make this discovery even more significant. 

The question that had to be answered concerning the axisymmetric 
instability was whether it was caused by the multiplicity of cells in the core 
or as a natural manifestation of the solid pipe wall boundary condition which 
is absent in the case of a trailing line vortex in an infinite domain. 
Figures 5.21a through 5.21d represent the eigenfunctions associated with this 
instability. For the sake of space, we will show only the amplitudes in this 

section. In figure 5.21a, it is shown that wall presence has resulted in a 

second peak in the radial component of the eigenfunction. A similar situation 
exists for the axial component of the eigenfunction except that the magnitude 
of the peak close to the wall is much less than that of the center peak. 
However, based on the shape of the eigenfunctions close to the centerline, the 
similarity between this mode and the axisymmetric viscous mode of Batchelor's 
vortex is striking. Overall this is a C type disturbance with most of the 
perturbation energy being concentrated near the centerline. The variation of 
with a for the axisymmetric disturbance is presented in figure 5.22. It 

must be mentioned that the effect of setting tp equal to zero on this mode was 

negligible. The variation of w. with Reynolds number for q = 1.0, and « = = 

2.0, is displayed in figure 5.23. The growth rate increases as Re increases, 
attaining its maximum value at moderate Reynolds numbers and then asymptoting 
to a constant value or dropping slightly. This mode clearly is an inviscid 
one. The variation of w r with Re is shown in figure 5.24. The curve is 
almost a straight line with constant slope. 

Figure 5.25 shows the variation of w. with a for an asymmetric mode (n = 
1) with q = 0.4. The first mode, which has much larger growth rates, is a 
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Amplitude 


5.21 Amplitude of eigenfunction components of axi symmentnc 
(n = 0) perturbation for a two cell vortex, a) radial 
velocity, b) azimuthal velocity, c) axial velocity, d) 

pressure. Here, a = 2.25, q = 1.0, Re = 10^, and 4 > ^ 0 
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Fig. 5.21b) Continued. 
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Fig. 5.21c) Continued. 
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Ref. 5.21d) Concluded. 
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short wave instability while the second mode exists only at the low end of the 
wavenumber spectrum (long wave instability). The effect of setting 4 = 0 on 
these modes is presented in figure 5.26. While there is hardly any change in 
the primary mode, note the significant drop in the growth rate of the second 
mode. The term associated with <\> represents the non-parallel effect in the 
velocity profiles of Donaldson and Sullivan [33]. The second mode is 
displayed in more detail in figure 5.27. The instability is clearly a long 
wave instability. Note that at this particular value of the Reynolds number 
(Re = 10000) the growth rates are orders of magnitude higher than the n - 1 
viscous mode associated with the trailing line vortex represented in figure 
4.25. The real part of frequency of the second mode is presented in figure 
5.28. The variation of ay with a is linear. It was found that ay did not 
vary with 4 . Setting 4 equal to zero, the eigenfunctions of the first and 
second modes were obtained. The eigenfunctions associated with the second 
mode (at a = 0.2) are displayed in figures 5.29a through 5.29d. The 

resemblance between these curves and the eigenfunctions for the asymmetric 
viscous mode of the trailing line vortex (figures 4.27a-4.27h) is striking. 
The only observable difference is the position of the peak in the axial 
component of the eigenfunction. Here the peak occurs close to r = 0.7, which 
is natural since, as mentioned earlier, the inflectional point in the axial 
component of mean velocity is situated very close to this value of r. 
However, the second mode is basically a C type perturbation. The eigenfunc- 
tions of the first mode (at a = 2.5) plotted in figures 5.30a through 5.30d 
revealed the existence of co-instability even though the growth rate curve 
does not show any multiple peaks and is smooth. At this value of a, the 

perturbation is a C R W type. Notice how different the structure of this 
disturbance is when compared to that of the second mode even though both have 
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c 07 Detailed representation of the variation of the growth rate 
5 - 27 of 1 thedsymmetric (n • . 1 . second mode) disturbances with 
wavenumber for a two cell vortex. Here, q - 0. , e 
and the influence of nonparallel effects (4- * 0) ls a,so 
show 
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Fiq. 5.29 Ajnplitude of eigenfunction components of asymmetric 

(n = + 1 , second mode) perturbation for a two cell vortex, 
a) radial velocity, b) azimuthal velocity, c) axial velocity, 

d) pressure. Here, a = 0.2, q = 0.4, Re = 10^, and 4 . = 0. 
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Fig. 5.29b) Continued. 
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Fig. 5.29c) Continued. 
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Fig. 5.29d) Concluded. 
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Fia. 5.30 Amplitude of eigenfunction components of asymmetric 
ng. d.ju nmp. + fipst perturbation for a two cell vortex. 

a) radial velocity, b) azimuthal velocity, c) axial velocity, 

d) pressure. Mere, a = 2.5, q = 0.4, Re = 10^. ancl 4 * = 0 
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the same value of n. The behavior of the growth rate at a higher value of q 

for the n = 1 disturbance is presented in figure 5.31. The primary mode 

indicates a slight damping effect. Figure 5.32 shows how ay of the primary 

mode varies with the axial wavenumber for q = 0.8. The inviscid nature of the 

first mode at q = 0.8 and a = 2.5 is shown in figure 5.33. The growth 

rate, w i , drops slightly as Re increases. However, calculations at much 

higher values of Re are needed to determine the asymptotic value of the growth 

rate, u. . The variation of the real part of frequency with Reynolds number 

is presented in figure 5.34. There is considerable change in the magnitude of 

a) over intermediate values of Re. The behavior of the second mode at this 
r 

higher value of swirl parameter (q = 0.8) is displayed in figure 5.35. The 

higher amount of swirl has resulted in higher growth rates as well as 

increasing the range of axial wavenumber, a, where instability exists. 

However, the instability still is a long wave one. The real part of frequency 

of the second mode varies very rapidly with a as shown in figure 5.36. The 

curve obtained with <\> = 0 is identical to figure 5.36 where <6^0. The 

viscous nature of the second mode is displayed in figure 5.37 for a = 0.5 and 

q = 0.8. Note that as the Reynolds number becomes large the effect of the <j> 

term is diminished. Also at high Re, the destabilizing effect of viscosity is 

greatly diminished. Figure 5.38 indicates little influence of Reynolds number 

on a) for the second mode, 
r 

The variation of the growth rate with a at q = 0.8 for the n = 2 mode is 
shown in figure 5.39. This is a long wave instability very much like the n = 
1 case presented in figure 5.27. For this mode, no sign of instability at 
high values of a was observed. The variation of ay with a for this mode is 
presented in figure 5.40. The effect of viscosity on the n = 2 mode at a = 
0.5 and q = 0.8 is displayed in figure 5.41. At high Reynolds numbers, the 
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Fig. 5.31 Variation of the growth rate of 
asymmetric (n = +1) disturbances 
with wavenumber for a two cell 
vortex. Here, q = 0.8, Re = 10^ 
and ii; * 0. 















Fig. 5.39 Influence of non-parallel effects 
on the variation of the growth 
rate of n = +2 disturbances with 
wavenumber for a two cell 
vortex. Here, q = 0.8, Re = 10^, 







ontr ibut ion from the 4 * term is small. However, the behavior at low Re is 
uite different. It was found that as the Reynolds number decreases, the 
urve associated with 4 , * 0 dipped to a minimum value. Even though not shown, 
urther decreases in the value of Re caused the growth rate, u> r to attain 
arge values, becoming infinite as Re + 0. The effect which was pointed out 
y Lessen (private communication) is an artifact of the non-parallel mean flow 
,nd means that the viscosity is causing the vortex core to spread at an 

infinite rate. By inspection of figure 5.41, the growth rate is asymptoting 
;o a constant value as Re increases which makes this mode more likely to be an 
inviscid one. This instability needs further testing. The dependency of u> p 
an Reynolds number for the n = 2 mode is displayed in figure 5.42. There is a 
slight increase in the value of the real part of the frequency as Re 

increases. The effect of setting c|> - 0 on was negligible and hence is 

not presented here. 

The variation of the growth rate with axial wavenumber at higher values 
of the swirl parameter (q = 1.2) is shown in figure 5.43. The added rotation 
has a destabilizing effect especially on the curve associated with 4 , = 0. The 
peak growth rate is slightly shifted to higher a with a net increase in the 
range of axial wavenumber where instabilities are present. Variation of w p 
with a at this q, shown in figure 5.44, indicates a similar trend to that of 
the q = 0.8 case. However, note that the slope of the present curve is much 
smaller in comparison to the case displayed in figure 5.40. 

The variation of growth rate, with a for the n = -1 mode at q - 0.4 
is shown in figure 5.45. The complexity of each individual curve shows signs 
of co-instabilities especially at the low end of the axial wavenumber 

spectrum. The maximum growth rates are twice as large as the growth rates 

associated with the trailing line vortex presented in figure 4.4. Also note 
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Fig. 5.43 Variation of the growth rate of 
n = +2 disturbances with 
wavenumber for a two cell vortex 
Here, q = 1.2, Re = 10 . 
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Fig- 5.44 Variation of the real part of 

frequency of n = +2 disturbances 
for a two cell vortex. Here, 
q = 1.2, Re = 10^, and ^ * 0. 
Identical curve is obtained when 

ij/ = 0 . 
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5.45 Variation of the growth rate of 
asymmetric (n = -1) disturbances 
with wavenumber for a two cell 
vortex. Here, q = 0.4, Re = 10 
and ijj * 0 . 



the enormously widened ranges of a where instabilities are present. It must 
be mentioned that higher modes were present in this case but they have been 
omitted for clarity. However, the most important feature of the present 
figure is the crossing of the curves associated with the first and second 
modes. This is the first time that such an effect has been observed. In all 
of the previous works (including Batchelor's vortex. Long's vortex and 
rotating Poiseuille pipe flow) each successively higher mode was entirely 
engulfed by the previous mode (e.g. figure 4.4). In such a flow, as has been 
pointed out by Maslowe and Stewartson [20], once the value of a exceeds the 
neutral mode value for the primary mode, the flow becomes stable to all such 
inv i sc id disturbances. This is not the situation here and it seems that this 
is a property of the multiple cell vortices. But the question that remains is 
why Foster and Duck [23] did not observe such an effect in their inviscid 
study of a two cell Long's vortex. Of course, it must be emphasized that 
their calculations were for a vortex in an infinite domain. 

Another point to be made here is that the crossings of the instability 
curves has made numbering of these modes quite arbitrary. Such difficulty has 
been encountered previously in the case of rotating pipe flows by Cotton and 
Salwen [18]. In the present section, we have numbered each mode according to 
their maximum growth rates. The variation of a> p with axial wavenumber for 
the first two modes are shown in figure 5.46. It seems that there is a 
degenerate point occurring near a * 0.8. At this particular value of a, the 
complex frequency for the first mode, u^, is approximately equal 

Although figure 5.46 does not indicate that such a state is reached exactly, 
many such degenerate points were encountered by Cotton and Salwen [18] in 
their viscous study of rotating Poiseuille pipe flow. The inviscid nature of 
the n = -1 mode at a = 0.6 and q = 0.4 is presented in figure 5.47. As Re 
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Fig. 5.46 Variation of the real part of 

frequency of asymmetric (n = -1) 
disturbances wi th wavenumber for a 
two cell vortex. Here, q = 0.4, 

Re = 10^, and ip * 0 . 
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Fiy. 5.47 Influence of Reynolds number on 
the growth rate of asymmetric 
(n = -1) disturbances for a two 
cell vortex. Here, a = 0.6, 
q = 0.4, and ij> t 0. 



increased beyond 10,000, the curves show little variation with the Reynolds 
number. Figure 5.48 shows a similar trend in the case of the real part of 
frequency, uy. The inviscid behavior of this mode at higher values 
of a (a = 2.5) is shown in figure 5.49. The variation of the growth rate 
with axial wavenumber for the n = -1 mode at q = 0.8 is shown in figure 
5.50. Again, the evidence of co-instabilities is very clear. The addition of 
swirl has had a strong destabilizing effect. However, one important 
feature is the resemblance of the instability curves for 0 < a < 1.0 to that 

of the trailing line vortex shown in figure 4.4. 

At this point, it is interesting to look at the nature of the co- 
instabilities (compound instability modes) associated with the primary mode. 
This task can be accomplished by moving in axial wavenumber space from high 
values of a towards the low end of the spectrum while looking at the shape of 
the eigenfunctions at selected values of a. For the sake of space, we are 
going to show only the radial component of the eigenfunction since this is the 
component most frequently dealt with in the open literature. It was found 
that the nature of radial eigenfunction is typical enough in highlighting the 
basic mechanism involved in co-instabilities so that we can omit displaying 
the other components of the disturbance eigenfunctions. Figure 5.51a displays 
the variation of the radial perturbation with radial distance at a - 4.0. 
Most of the energy of the wave is concentrated at the wall and center part of 
the flow with the peak occurring at r ~ 0.6. For this axial wavenumber, 
the disturbance is an R W type. At a = 3.0, as shown in figure 5.51b, a 

small peak in the eigenfunction has appeared at r = 0. However, most of the 
energy remains in the outer part of the flow and still the perturbation can be 
considered an R W type. But the situation is quite different at o = 2.6 which 
is shown in figure 5.51c. Here, the center mode has become the dominant mode 
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Fiy. 5.4b Influence of Reynolds number on 
the real part of frequency of 
asymmetric (n = -1) disturbances 
for a two cell vortex. Here, 
a = 0.6, q = 0.4, and ib * 0 . 




49 Influence of Reynolds number on 
the growth rate of asymmetric 
(n = -1) disturbances for a two 
cell vortex. Mere, a = 2.5, 
n = D.4. and il> *■ 0 . 



Fig. 5.50 Variation of the growth rate of 
asymmetric (n = -1) disturbances 
with wavenumber for a two cell 
vortex. Here, q = 0.8, Re = 10 4 
and <l> * 0 
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Fin 5 51 AmDlitude of radial eigenfunction component of an asymmetric 
Flg * ' ( n p =_i > first mode) perturbation for a two cell vortex. 

a) a = 4.0, b) ot = 3.0, c) a = 2.6, d) a = 2.0, 

e) a = 1.1, f) a = 0.4. Here, q = 0.8. Re = 10 , and * 0 . 
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although most of the energy now is distributed over the entire pipe radius. 
This disturbance is clearly aCRW type. At a = 2.0 (shown in figure 

5.51d), the peaks associated with R and U have subsided and C has become the 
only significant instability zone. The variation of the eigenfunction with 
radial distance for a = 1.1 is presented in figure 5.51e. At this axial 
wavenumber, the two peaks in the outer flow part have almost entirely dis- 
appeared. At a = 0.4 (presented in figure 5.51f), the C type instability 
is the sole mode present with no sign of R or W type modes. For this 

value of a, the radial eigenfunction approximates the shape of the radial 
component of the n = 1, long wave instability (figure 4.27a) of a Batchelor 
vortex. 

The variation of growth rate, u., with axial wavenumber, a , at q = 0.8 
for the n = -2 mode is shown in figure 5.52. There are signs of co-instabili- 
ties for each curve. For this disturbance, the range of a where instabilities 
are present, has doubled. But notice the huge growth rates of the primary 
mode. The maximum value of for the first mode indicates that the 

disturbance has very short e-folding time. Variation of the real part of 
frequency with a for the first two modes is shown in figure 5.53. As a + 0 
the problem of mode jumping is very much apparent from figures 5.52 and 5.53. 

The behavior of n = -2 perturbations with Reynolds number at q = 0.8 and a = 

1.0 is shown in figure 5.54. The growth rate is independent of Re for 
Reynolds numbers as low as 4000. 

It must be clear by now that the behavior of most disturbances are 
generally different at low values of a than those associated with high axial 
wavenumbers. It seems that almost all instability modes termed as viscous 

modes are only present in the form of long wave instabilities; while evidently 

the inviscid modes are only manifested in the form of intermediate or short 
wave instabil ities. 


174 



Amplitude 


1.0 



Fig. 5 . 5 Id ) Continued. 
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Fig. 5.51e) Continued. 
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Fig. 5.51f) Concluded. 
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Fig. 5.52 Variation of the growth rate of 
n = -2 disturbances with 
wavenumber for a two cell 
vortex. Here, q = 0.8, Re = 10 
and ^ * 0. 
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Fig. 5.53 Variation of the real part of 

frequency of n s -2 disturbances 
with wavenumber for a two cell 
vortex. Here, q = 0.8, Re = 10 , 
and tjf ^ 0 . 
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5.3 Stability of Three Cell Vortex at Re = 4000 


Due to space limitations in this section only the variation of growth 
rates with a for different values of azimuthal wavenumbers are presented. In 
the case of the n = -4 mode, the effect of increased swirl on the instability 
is shown clearly. However, a detailed presentation for a three cell vortex 
similar in extent to that of Section 5.2 is deferred to future studies. 

The three velocity components for a three cell vortex are presented in 
figure 5.55. Note the core expansion in the azimuthal velocity component 
which is caused by the flow reversal of the axial velocity on the 
centerline. This flow reversal has displaced the vortex from the central 
region and depleted the core area of any rotational kinetic energy. Hence, we 
expect non-parallel effects to play an important role, especially at finite 
Reynolds numbers. 

The variation of the growth rate, u. , with <* for the axisymmetric mode, n 
= 0, at q = 0.19 is shown in figure 5.56. The first mode is very similar in 
extent to the axisymmetric disturbance for a two cell vortex as displayed in 
figure 5.20. However, the second mode is peculiar since the growth rate 
asymptotes to a constant as a approaches zero. Further testing revealed that 
the second mode of instability appeared to be caused by the same mechanism 
discussed in the previous section in relation to figure 5.41 To test this 
observation, the growth rate curve was obtained again while setting <|> = 0, 
as shown in figure 5.57. The growth rate of the second mode became so small 
that it could not be plotted on the same scale as the first mode. It must be 
emphasized here that further testing is needed in order to understand the 
nature of this second mode of instability. Of special importance is the 
variation of the growth rate with Reynolds number. 
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Fig. 5.55 Velocity components for a three 
cell vortex. 
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5.56 Variation of the growth rate of 

axi symmetri c (n = 0) disturbances 
with wavenumber for a three cell 
vortex. Here, q = 0.19, 
lie = 4000, and # 0* 





The variation of the growth rate with axial wavenumber for the asymmetric 
mode, n = 1, at q = 0.19 is shown in figure 5.58. The primary mode has 
enormous growth rates over a wide range of a. Also there is evidence of co- 
instabilities for the second mode. The nature of the third mode (long wave 
instability) makes it a prime candidate to be a viscous mode. The non- 
parallel term (<p) was set equal to zero and the results are shown in figure 
5.59. It is important to emphasize here that mode three did not disappear 
altogether, rather its amplification rate, u . , became very small, preventing 
its plot on the same scale as the other two modes. However, as in the case of 
axisymmetric mode, further testing is required to understand the nature of 
this instability. The influence of on the first and second modes is 
negligible. Actually, in the case of n < 0, the parameter <j> has no drastic 
effect on the nature of instabilities and hence its influence will not be 
discussed any further in this section. 

Figure 5.60 displays the variation of the growth rate with a for the n 
= -1 mode at q = 0.19. For this value of azimuthal wavenumber, there are two 
modes with comparable growth rates. Notice how different these curves are 
when compared to the curves obtained for the two cell vortex. Incidentally, 
the eigenfunctions revealed that these modes are located near the two 
inflectional points present in the axial velocity profile. 

The behavior of growth rate with a for an n = -2 disturbance at q - 0.19 
is shown in figure 5.61. The primary mode has a growth rate comparable to 
that of n = -1 perturbations as presented in the previous figure. There are 
signs of co-instabilities especially for the first mode. The variation of ^ 
with axial wavenumber for the n = -3 disturbance at a swirl number of 0.19 is 
displayed in figure 5.62. The existence of co-instabilities is very evident. 
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Fig. 5.58 Variation of the growth rate of 
asymmetric (n = +1) disturbances 
with wavenumber for a three cell 
vortex. Here, q = 0.19, 

Re = 4000, and d> * 0. 
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Fiq. 5.59 Variation of the growth rate of 
asymmetric (n = +1) disturbance 
with wavenumber for a three cell 
vortex. Here, q = 0.19, 

Re = 4000, and iji = 0. 
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Fig. 5.61 Variation of the growth rate of 
n = -2 disturbances with 
wavenumber for a three cell 
vortex. Here, q = 0.19, 

Re = 4000, and * 0. 
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Figure 5.63 displays the variation of the growth rate with a for an n - 
-4 perturbation at q = 0.19. Note how similar the behavior of the instability 
curves are compared to those of the n = -3 mode shown in the previous 
figure. Again, the co-instabilities are present for all three curves plotted 
in the figure. Figure 5.64 presents the results obtained at q = 0.38. The 
increase in swirl has a strong destabilizing influence on all of the modes. 
It has also increased the effective range of a where instabilities are 
present; and crossing of the second and third modes is noted. Further 
increases in the swirl parameter to a value of 0.57 (shown in figure 5.65) has 
resulted in even larger growth rates. The crossing of the second and third 
modes is pronounced. Additional increases in the swirl number to a value of 
0.76 produces further increases in the growth rates. The results which are 
shown in figure 5.66 indicate that the primary mode's growth rate has very 
short e-folding time, attaining a maximum value of 1.22 at a = 4.8. Also, 
preliminary calculations indicate that in the cases of multiple cell vortices 
substantially more swirl is needed to stablize the flow with respect to these 
inviscid disturbances than is needed for a trailing line vortex. 

The variation of the growth rate with axial wavenumber at q = 0.19 for an 
n = -5 perturbation is shown in figure 5.67. The curves show very similar 
behavior to those modes with lower negative azimthual wavenumbers. Finally, 
the variation of the growth rate with a for an n = -10 mode at q - 0.19 is 
presented in figure 5.68. The range of a over which the primary mode 
remains unstable has increased considerably. This is very much in line with 
the asymptotic analysis of [20, 26, 29-31] which predict that at sufficiently 
large negative n, the proper scaling for the disturbance axial wavenumber 
becomes a/n. 
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Fiy. 5.63 Variation of the growth rate of 
n = -4 disturbances with 
wavenumber for a three cell 
vortex. Here, q = 0.19, 

Re = 4000, and ii> * 0. 
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Fig. 5.64 Variation of the growth rate of 
n = -4 disturbances with 
wavenumber for a three cell 
vortex. Here, q = 0.38, 

Re = 4000, and ^ it 0. 
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Fig. 5.65 Variation of the growth rate of 
n = -4 disturbances with 
wavenumber for a three cell 
vortex. Here, q = 0.57, 

Re = 4000, and i|< * 0. 
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Fig. 5.66 Variation of the growth rate of 
n = -4 disturbances with 
wavenumber for a three cell 
vortex. Here, q = 0.76, 

Re = 4000, and tj> * 0. 
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Fig. 5.68 Variation of the growth rate of 
n = -10 disturbances with 
wavenumber for a three cell 
vortex. Here, q = 0.19, 

Re = 4000, and p * 0. 



One interesting observable fact here is the existence of so many 
significant unstable disturbances for different values of n at such a low 
value of the swirl parameter (q = 0.19). Also, the relatively equal strength 
of the various primary instabilities (we are speaking of the primary modes for 
different negative n) at this constant q is a bit surprising too. Both of the 
above features are not present for a Batchelor's vortex. 


198 



Chapter 6 


CONCLUSIONS 

The present study was undertaken to obtain the effects of initial mean 
velocity profile on the stability characteristics of a streamwise vortex with 
the objective of developing control and manipulation strategies for 
engineering applications. To that end, the temporal stability of multiple 
cell vortices was considered. It was shown that cell multiplicity in the 
vortex core has drastic effects on the stability characteristics. That is the 
multiple cell structure is responsive to a wide range of disturbances with 
different azimuthal wavenumbers. This effect is most pronounced in the case 
of n > 0 modes where previous studies of single cell trailing line vortices 
(as well as the two cell Long vortex) had indicated stability. Here, it was 
shown that multiple cell vortices are highly unstable, at least for n = 0, +1, 
+2 modes, in addition to the usual negative n modes. The enormous growth 
rates obtained for the axisymmetric mode (n = 0) was surprising since this 
particular type of disturbance has never been observed to be unstable before, 
even in the case of a two cell Long vortex. Overall, on the basis of the 
numerical calculations, it is concluded that the growth rates of instabilities 
in multiple cell vortices are substantially larger (two to three fold 
increases are observed) than those of a single cell trailing line vortex. It 
was also found that there is a substantial increase in the effective range of 
axial wavenumbers where instabilities are present in multiple cell vortices. 
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Another important feature whose existence has been shown in the present 
work is co-instability (short for compound instability modes). These 
instability modes have eigenfunctions which show features that can be 
characterized as center- and ring-modes as well as wall mode instability zones 
simultaneously. In the case of two and three cell vortices, the growth rate 
curves (for n < 0 only) possess multiple peaks (lobes) because of the co- 
instabilities. Although the existence of co-instability has been conjectured 
by previous asymptotic analyses [20, 26, 32], they had never been produced 
numerically. 

In addition, this work has shown that viscous forces play a far greater 
role in vortex stability than had been believed previously. The variety of 
viscous modes obtained for the Batchelor vortex as well as for multiple cell 
vortices showed that any future stability calculations must include the 
effects of viscosity. Although this subtle role of viscosity was suspected by 
Masl owe and Stewartson [20], its extent and influence were unknown prior to 
this study. 

These viscous disturbances, without any exceptions are long wave 
instabilities. Their wavelengths scale with the core radius and have maximum 
growth rates which are orders of magnitude smaller than the inviscid modes. 
However, based on experimental observation and measurements, these viscous 
instability modes appear to be the only physically observed modes. 

Based on the present numerical calculations, multiple cell vortices have 
desirable instability charater istics in that breakup and dissipation of these 
types of trailing line vortices can be accelerated. Clearly, it is not 
necessary to alter the radial variation of the azimuthal component of velocity 
and the radial component of the velocity is quite small and probably would not 


200 



be controlled. Only the axial velocity component must be modified and it is 
necessary to induce at least a two cell structure in the flow. The desired 
flow is downstream jet along the vortex axis, surrounded by an upstream flow 
in the outer portion of the vortex. This opportunity for accelerating vortex 
instability warrants further study. 

The Chebyshev spectral collocation algorithm, with staggered grid, was 
developed to study the stability of multiple cell vortices. The method is an 
extension of the technique developed by Khorrami et al . [10]. The numerical 
stability problem has been formulated in primitive variable form, evaluating 
the velocity components at the grid points and staggering the pressure at the 
mid-grid points. The staggered pressure approach has eliminated the need for 
artificial pressure boundary conditions. Direct comparison with other 
formulations has shown that the method is robust and produces accurate 
results. However, comparison with the non-staggered spectral method of 
Khorrami et al . [10] has shown that contrary to prevailing opinion, when the 
artificial pressure boundary conditions are employed correctly they do not 
affect the accuracy or the convergence rate of the spectral methods. 

While validating the accuracy and convergence behavior of the present 
method, two new, long wavelength viscous modes of instability for Batchelor's 
vortex were found. These two forms of instabilities include an axisymmetric , 
n = 0, as well as an asymmetric, n = 1 mode. This work has shown that the 
growth rates associated with these viscous disturbances are orders of 
magnitude less than the explosive growth rates of the inviscid perturba- 
tions. However, examination of contrail photographs [58, 62] has shown that 
both of these new modes can be identified in actual large trailing line 
vortices as well as in experimental measurements conducted in laboratory wind 
tunnels, [60, 61, 62]. Based on the present results, the following 
conclusions can be drawn: 
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1. Batchelor's vortex represents a good approximation to trailing line 
vortices and can be employed for further stability studies. 

2. Although the stability calculations were conducted for a laminar 

vortex, the excellent agreement obtained with experimental 

measurements and contrail photographs, which involve mostly turbulent 
vortices, suggests that these long wave viscous instability modes, 

once activated, do not interact significantly with the surrounding 
turbul ent field. 

Finally, the numerical method used in this study can be adopted easily to 
include cases where density stratification or compressibility effects play an 
important role. From the enormous body of data on vortices in turbulent 
atmospheres, the role played by stratification on vortex behavior is 

significant. Since trailing line vortices coexist in such a background, the 
effects of density stratification on the instability modes must be examined. 
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